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1  Overview  of  the  Research  Effort 


The  objective  of  the  current  research  is  to  develop  a  probabilistic  analysis  framework  incorporating  high 
performance  composite  design  methodologies  to  predict  the  reliability  of  composite  structures.  The  research 
links  state-of-the-art  laminate  design  techniques  with  failure  models.  This  effort  has  the  potential  for 
application  to  the  types  of  engineering  tasks  that  are  required  for  the  development  of  new  rotorcraft  using 
modem  composite  materials  and  technologies  and  is  especially  relevant  to  the  areas  of  damage  tolerance  and 
crash-worthiness.  The  research  will  yield  a  methodology  and  software  that  would  greatly  improve  both  our 
physical  understanding  and  analytical  capability  in  this  subject  area. 

The  methodology  is  being  demonstrated  on  a  helicopter  rotor  hub  test  specimen  with  thick,  thin,  and 
tapered  regions.  The  specimen  is  subjected  to  a  static  axial  load  P  and  an  oscillatory  angular  displacement  <l>. 
The  failure  mechanism  that  is  observed  from  this  type  of  fatigue  loading  is  an  initial  tension  crack  followed 
by  internal  delamination  at  the  thick  to  taper  transition  where  internal  ply  drop-offs  occur. 

A  stmctural  model  of  the  rotor  hub  was  created  using  the  finite  element  program  ANSYS.  The  output 
from  the  stmctural  analysis  in  ANSYS  was  fed  into  the  Virtual  Crack  Closure  Technique  (VCCT)  module  to 
model  del^ination  onset.  A  Reliability  Analysis  framework  including  Response  Surface  Analysis  and 
FORM  modules  was  developed  to  integrate  with  the  stmctural  analysis  modules.  The  response  surface 
analysis  was  used  to  develop  a  limit  state  equation  relating  the  primitive  input  design  parameters  (random 
variables  such  as  P,  and  <!))  to  strain  energy  release  rate  FORM  analysis  estimated  the  probability  of 
delamination  onset  and  furthermore,  determined  the  distribution  of  the  failure  probability  over  a  range  of 
cyclic  lives.  In  addition,  the  sensitivity  analysis  in  FORM  revealed  that  the  external  loading  to  the  stmcture, 
namely  P  and  O  are  the  most  critical  parameters  affecting  the  reliability  of  the  rotor  hub.  Summarizing,  the 
following  tasks  were  successfully  completed  at  the  end  of  Phase  1. 

■  Demonstrated  probabilistic  analysis  framework 

■  Demonstrated  composite  Finite  Element  modeling  (FEM)  technique 

■  Developed  framework  for  integrating  probabilistic  and  FEM  codes 

■  Applied  developed  method  to  a  practical  problem 

Phase  I  of  the  proposed  research  effort  has  successfully  demonstrated  that  the  probabilistic  analysis 
framework  can  be  integrated  with  current  state-of-the-art  composite  design  to  estimate  the  failure  probability 
due  to  various  critical  failure  modes.  In  addition,  the  same  framework  can  be  used  to  compute  the  sensitivity 
of  the  reliability  of  the  final  design  to  the  random  variables  without  any  extra  computational  effort. 


2  Identiffcation  and  Significance  of  the  Problem 

The  structural  evaluation  of  composite  materials  is  different  from  that  used  in  isotropic  materials. 
Behavior  of  composites  is  dependent  on  the  specific  way  in  which  a  given  composite  is  layered.  The 
orientations  allow  the  designer  freedom  in  tailoring  structures  for  obtaining  the  desired  response.  The 
evaluation  of  composite  structures  needs  to  account  for  the  effects  of  tailoring. 

The  development  of  advanced  high-performance  composites,  until  now,  has  mainly  focused  in  achieving 
high, modulus  and  high  strength.  But  along  with  the  high  strength,  such  materials  must  also  be  able  to  absorb 
energy  and  resist  cyclic  fatigue  loads.  Fatigue  loads  can  initiate  progressive  failure  in  a  composite  laminate  by 
way  of  successive  delamination,  matrix  cracking,  fiber  waviness,  etc.  Such  imperfections  can  severely  reduce 
the  load  canying  capacity  of  the  laminates,  ultimately  leading  to  structural  failure.  While  it  is  possible  to 
address  the  effect  of  geometric  and  material  non-linearity  using  current  sophisticated  analysis  methods,  the 
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highly  stochastic  nature  of  damage  progression  and  extreme  sensitivity  to  small  geometrical  uncertainties 
results  in  untrustworthy  analytical  predictions. 

The  problem  is  further  compounded  by  the  fact  that  there  is  significant  randomness  in  the  material  and 
geometric  properties  of  such  composites.  Another  vital  source  of  variance  is  the  loads  to  which  composite 
structures  are  subjected.  The  exact  load  distribution,  its  variation  in  time  and  therefore,  its  structural  effect  on 
composites  for  aerostructure  applications  are  seldom  deterministic.  A  probabilistic  analysis  framework 
incorporating  the  current  laminate  design  methodologies  is  the  best  strategy  to  account  for  these  variations  in 
both  load  and  material  properties. 

3  Technical  Objectives  of  the  Phase  I  Effort 

The  research  developed  a  probabilistic  analysis  framework  incorporating  current  aerostructure  industry 
composite  laminate  design  methodologies  to  predict  the  reliability  of  composite  structures.  The  research 
considered  the  reliability-based  composite  design  methodology  as  being  divided  into  two  distinct  but  closely 
coupled  modeling  techniques.  The  fu^t  modeling  technique  is  the  structural  model  that  uses  finite  element 
analysis  to  determine  the  global  and  ply  level  response  of  the  structure  at  the  damage  sites.  The  finite  element 
analysis  is  the  same  as  that  currently  employed  in  the  aerostructure  industry.  The  finite  element  model  is 
extended  to  address  the  random  nature  of  the  structural  application.  The  model  considers  the  stochastic 
aspects  of  the  input  loads.  Variation  in  the  global  geometiy  is  considered,  along  with  local  ply  level  variation. 
The  distributions  of  geometric  variations  were  determined  by  investigating  the  sources  of  the  variations  such 
as  standard  manufacturing  techniques,  manufacturing  anomalies,  and  service  duty. 

The  second  modeling  techniques  are  the  failure  models  that  are  closely  coupled  with  the  structural  model. 
The  failure  models,  which  address  initial  damage,  are  specific  models  that  have  been  developed  by  the 
rotorcraft  industry  and  the  Army  Research  Laboratory.  The  failure  models  consider  many  of  the  same  input 
variations  as  the  structural  models.  The  failure  models  also  consider  variations  in  the  local  responses  to 
determine  the  probability  of  occurrence  of  damage  states. 

The  objective  of  the  Phase  I  research  was  to  develop  and  validate  an  overall  progressive  failure 
methodology  using  the  approach  outlined  above. 

4  Research  Team 

The  research  team  comprises  of  PerSyst  Development,  Vanderbilt  University  and  Bell  Helicopter.  For 
the  probabilistic  methodologies  to  have  commercial  application,  these  methodologies  must  be  developed  with 
a  thorough  understanding  of  current  rotorcraft  design  practices.  Thus,  a  close  working  relationship  has  been 
fostered  between  PerSystA^anderbilt  and  Bell  Helicopter.  Bell  has  identified  rotor  fatigue  loading  as  a 
significant  probabilistic  design  consideration  due  to  the  large  amount  of  scatter  in  the  observed  fatigue  life  of 
the  large  scale  composite  rotor  test  specimens.  The  Phase  I  effort  is  targeted  at  developing  and  applying  the 
probabilistic  methodology  to  the  design  and  analysis  of  the  rotor  test  specimens  currently  used  by  Bell.  The 
finite  element-based  macrostructural  model  is  consistent  with  current  modeling  methods  used  at  Bell  and  the 
microstructural  fatigue  damage  accumulation  model  is  based  on  the  models  developed  through  the 
cooperative  efforts  of  Bell  and  the  U.S.  Army  Research  Laboratoiy.  Bell  is  also  helping  the 
PerSystA^anderbilt  team  identify  the  work  that  must  be  completed  in  the  Phase  II  effort.  TTiis  collaboration 
will  assure  that  the  final  product  addresses  the  current  and  near  term  needs  of  composite  rotorcraft  design  and 
assure  that  the  fmal  product  can  be  successfully  integrated  into  the  design  practice. 

f 

5  Phase  1  Research  Conducted 

The  focus  of  the  research  has  been  in  three  areas;  1)  determining  the  failure  mechanism  and  developing 
the  structural  models  to  simulate  the  local  conditions  governing  the  failure  mechanism  for  the  demonstration 
problem,  2)  developing  the  reliability  computation  framework  for  the  demonstration  problem  and  3) 
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developing  additional  reliability  algorithms.  The  focus  areas  have  been  further  subdivided  into  several  tasks. 
Each  of  these  tasks  is  discussed  below. 

5.1  Focus  Area  1:  Failure  Mechanisms  and  Models 


Hingeless,  bearingless  helicopter  rotor  hubs  are  being  designed 
using  laminated  composite  materials  to  reduce  weight,  drag,  and  the 
number  of  parts  in  the  hub  as  shown  in  Figure  1.  An  effective  elastic 
hinge  is  designed  integrally  to  the  composite  rotor  yoke  by 
incorporating  a  tapered  region  between  thick  and  thin  regions.  The 
varying  thickness  of  the  tapered  region  is  accomplished  by  dropping 
internal  plies.  The  thick-taper-thin  geometry  is  tailored  to  give  the 
proper  flapping  flexure.  An  example  of  a  composite  rotor  yoke  is 
shown  in  Figure  2. 

During  flight,  the  rotor  hub  yoke  experiences  axial  tension  loads 
as  well  as  bending  loads.  The  axial  tensile  loads  are  centrifugal 
loads  caused  by  the  rotation  of  the  rotor.  The  bending  loads  are 
caused  by  the  interaction  of 
the  rotor  passage  with  the 
fuselage  or  other  static  structures  and  aerodynamic  gust  loads.  The 
actual  load  experience  in  service  is  a  complex  random  spectrum 
dependent  on  mission,  local  environment,  and  turbulence  [1]. 

Generalized  test  specimens  have  been  developed  to  understand 
the  basic  response  of  the  composite  rotor  hub  yokes.  The  specimens 
are  geometrically  simple  with  thick,  thin,  and  tapered  regions,  as 
shown  in  Figure  3  and  approximate  the  geometry  of  section  A-A  in 
Figure  2.  The  specimens  are  subjected  to  constant  axial  tensile  load 
(P)  to  simulate  the  rotation  of  the  rotor  at  constant  speed  and  cyclic 
bending  load  (V)  to  simulate  the  interaction  of  the  rotor  passage  with 
the  fuselage.  A  typical  test  setup  is  diagramed  in  Figure  3.  The  cyclic 
load  (V)  induces  an  angular  displacement  (6)  which  simulates  the 
flexural  bending  in  the  yoke. 


5. 7. 7  Demonstration  component 


Figure  1:  Rotor  hub  assembly. 
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5. 1. 2  Failttre  mechanism 

The  plydrop  in  the  laminate  creates 
geometrical  and  material  discontinuities 
that  cause  large  interlaminar  stresses  and 
initiate  delamination.  The  failure 
mechanism  observed  from  this  type  of 
loading  is  an  initial  tension  crack 
followed  by  internal  delamination  at  the 
thick  to  taper  transition  where  internal 
ply  drop-offs  occur. 

5.1.3  Failttre  model 

Several  studies  have  investigated  the 
effect  of  delamination  failure  in  tapered 
laminates.  Some  of  these  studies  have 
used  stress-based  criteria  for  modeling 
the  failure  [2-7].  Others  have  used  a 
strain  energy  release  rate  approach  to 
study  delamination  [8-13].  Most  of  the 
above  mentioned  methods  have  only 
considered  delamination  under  pure 
tension,  bending  or  torsion  loads.  Very 
few  studies  have  considered  the 
combined  effect  of  bending  and  tension 
on  delamination  of  tapered  laminated 
Figure  4:  Model  of  crack  tip  at  plydrop.  composites  [14,15]. 

In  this  study,  the  method  used  to 

model  the  delamination  failure  mechanism  is  based  on  total  strain  energy  release  rate  G  which  is  the  sum  of 
the  mode  I  strain  energy  release  rate  and  the  mode  II  stain  energy  release  rate.  A  virtual  crack  closure 
technique  (VCCT)  is  used  to  calculated  G  at  the  delamination  tip  as  shown  in  Figure  4  [16,17]  such  that 

G  =  Gf  +  Gji  (4.1) 

where, 

and  u  and  v  are  tangential  and  perpendicular  nodal  displacements  respectively  and  F,  and  F„  are  the  tangential 
and  perpendicular  nodal  forces  respectively 

Delamination  onset  is  assumed  to  occur  when  the  calculated  G  exceeds  the  G^^,  derived  from  material 
coupon  delamination  tests  [8,13,15].  The  details  of  estimating  the  probability  of  delamination  onset  are 
explained  in  the  Section  5.2.2.  A  finite  element  model  will  be  used  to  determine  the  local  forces  and 
displacements  needed  to  calculate  G  as  shown  in  Figure  4.  The  finite  element  model  is  discussed  below. 
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5.1.4  Structured  model 


The  geometry  of  the  ANSYS  finite  element  model  for  the  problem  is  shown  in  Figure  5.  The  model 
consists  of  6814  nodes  and  2860  elements.  The  elements  are  of  three  types:  4-node  quadrilateral,  8-node 
quadrilateral,  and  6-node  triangular  elements.  In  the  thick  to  taper  transition  region  where  delamination  is 
critical,  the  finite  element  mesh  is  refined  (see  Figure  6)  and  modeled  with  8-node  quadrilateral  elements.  The 
average  element  size  in  this  region  is  1/4*  the  ply  thickness.  Such  small  element  size  ensures  that 
computation  of  the  strain  energy  release  rate,  G,  does  not  suffer  from  numerical  inconsistencies.  The  resin 
pockets  at  the  end  of  the  dropped  ply  zones  are  modeled  with  6-node  triangular  elements.  The  tension  crack 
at  the  interface  of  the  dropped  ply  zone  and  the  resin  pocket  has  been  modeled  with  the  help  of  contact 
elements.  The  delamination  cracks  at  the  interface  of  the  belt  plies  and  the  dropped  ply  zones  have  been 
modeled  with  the  help  of  duplicate  multi  point  constraint  (MFC)  nodes. 

5.1.5  Calculation  of  laminate  material  properties 

A  simple  FORTRAN  code  has  been  written  to  compute  the  laminate  material  properties  fi’om  the  basic 
material  properties.  The  method  is  well  known,  and  is  described  in  any  textbook  on  composite  materials 
mechanics,  e.g.,  [18]. 

5.1.6  ANSYS  finite  element  stress  analysis 

A  nonlinear  finite  element  analysis  is  executed  for  the  stress  analysis.  The  total  load  increment  is  divided 
into  ten  steps.  In  each  step,  a  tenth  of  total  loading  is  applied.  The  sti^ess  matrix  is  reassembled  at  each  step 
and  a  maximum  of  30  iterations  is  used  for  convergence  of  the  stress  analysis.  A  batch  file  in  ANSYS  is 
written  to  extract  the  displacements  and  the  nodal  forces  at  the  elements  around  the  delamination  tip  as  a 
result  of  the  stress  analysis. 

A  program  in  C  language  is  written  to  use  the  displacements  and  the  nodal  forces  at  the  delamination  tip 
to  compute  the  strain  energy  release  rate,  G,  according  to  the  Eq.  (4.1).  The  various  steps  in  ANSYS  stress 
analysis,  including  the  building  of  the  response  surface  model  are  illustrated  in  a  flow  chart  in  Figure  7. 


in 


Figure  5:  Model  of  tapered  rotor  hub 
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Figure  6  Refined  mesh  in  the  thick  to  taper  transition. 
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5. 1. 7  Delamination  growth  analysis 

The  objective  of  this  task  is  to  calculate  the  strain  energy  release  rate  for  different  delamination  lengths 
and  determine  the  delamination  length  at  which  the  maximum  strain  energy  release  rate  {Gmad  is  observed. 
Then  the  response  surface  method  and  direct  FORM  can  be  used  to  calculate  the  delamination  probability  of 
helicopter  rotor. 

Figure  6  shows  the  detail  of  delamination  growth  analysis:  From  point  A  to  B,  the  duplicate  nodes 
between  elements  along  the  delamination  front  are  released  step  by  step.  In  each  step,  nonlinear  finite  element 
analysis  is  performed  and  the  G  value  at  the  delamination  tip  is  calculated.  Table  1  and  Figure  8  show  the 
relationship  between  the  delamination  length  and  the  corresponding  strain  energy  release  rate.  (Note: 
Delamination  starts  from  X  =  5.0,  and  grows  to  the  left).  The  maximum  G  value  of  1 .014553  in-lb/in^  is 
observed  when  the  delamination  length  equals  0.026167  inch.  The  model  with  this  delamination  length  is 
used  ftuther  for  response  surface  analysis  as  well  as  direct  FORM  analysis. 

Table  1:  G  value  and  delamination  length. 


Delamination 
length  (in) 

0.016033 

0.0211 

0.026167 

0.031234 

0.036302 

0.074312 

G  value 
(in-lb/in^) 

0.605323 

0.891599 

1.014553 

0.944279 

0.87124 

0.447136 

Delamination 
length  (in) 

0.114855 

0.195935 

0.25675 

0.317564 

0.378375 

4.581083 

G  value 
(in-lb/in^) 

0.258372 

0.17699 

0.192102 

0.237581 

0.336777 

0.459771 
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Figure  8;  G  values  for  different  delamination  lengths. 

5. 1.8  Random  variable  definition 

Based  on  data  from  Bell  Helicopter,  seven  random  variables  have  been  identified  for  probabilistic 
analysis  with  the  finite  element  model  of  the  helicopter  rotor  hub.  These  are  material  property  variables  En, 
E22,  v/j,  and  G/j,  the  oscillatory  bending  angle  6  caused  by  the  cyclic  loading,  the  magnitude  of  the  axial  load 
P,  and  the  critical  value  of  the  strain  energy  release  rate  Gcru.  All  the  random  variables  are  assumed  to  be 
Gaussian  (normal)  for  the  purpose  of  demonstration,  and  their  mean  values  and  standard  deviations  are 
shown  in  Table  2. 


Table  2  :  Random  variable  description 


S2/E7T1-2  Tape 

Property 

Mean 

Std.  Deviation 

Eu  (Msi) 

6.90 

0.09 

E22  (Msi) 

1.83 

0.05 

Gn  (Msi) 

0.698 

0.015 

0.28 

0.01 

P(kips) 

30.8 

3.08 

0  (degrees) 

12 

1.67 
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The  objective  of  the  demonstration  problem  is  to  estimate  the  probability  of  initiation  of  delamination  at 
the  required  life.  This  is  assumed  to  occur  when  the  strain  energy  release  rate  exceeds  the  limiting  value  Gent, 
which  is  a  function  of  load  cycles,  N.  Figure  9  shows  the  Gcru  vs.  as  supplied  by  Bell  Helicopter  based  on 
test  data.  Figure  9  also  shows  the  best  fit  line  through  the  data  assuming  Gent  is  a  linear  function  of  Logio(N). 
According  to  Figure  9,  the  relationship  between  Gc„v  and  iV  can  be  expressed  as: 

G„,,  =  448.56  -  58.57  liog,o(iV)  (4.2) 

Based  on  statistical  analysis  of  the  Gcrn  vs.  A^data,  it  was  determined  that  Gent  is  a  Gaussian  (normal) 
random  variable  with  its  mean  value  described  by  Eq.  (4.2)  and  a  standard  deviation  of  36.6  J/m^. 


♦  Data 

- Mean 

. +  Sigma 

- Sigma 


Figure  9:  Relationship  between  G_crit  and  cycles  of  life. 


5.2  Focus  Area  2:  Development  of  Reliability  Analysis  Module 

This  task  requires  developing  a  reliability  analysis  module  using  the  First  Order  Reliability  Method 
(FORM)  [19].  The  methodology  of  FORM  requires  describing  the  failure  condition  in  a  mathematical  form 
as  follows: 

S  =  Gcrii  -  (4.3) 

Eq.  (4.3)  is  called  the  limit  state  function.  The  limit  state  is  a  characteristic  of  the  system  component  that 
can  be  observed,  and  is  defined  to  represent  the  failure  of  the  component  to  attain  a  required  level  of 
performance.  According  to  Eq.  (4.3),  g  ^  0 ,  i.e.,  G»ra»  exceeds  GcHh  represents  failure  and  g  >  0  represents 
the  safe  state  of  the  system.  The  limit  state  function  relates  the  system  performance  to  primitive  random 
design  variables.  In  this  case,  Gan  and  G„atare  random  variables  that  are  dependent  on  primitive  random 
variables  such  as  material  properties,  geometric  properties,  loading  conditions,  etc.  which  have  inherent 
scatter  in  their  definition.  Gaii  can  be  considered  to  be  the  capacity  of  the  system  and  Gmarcan  be  considered 
to  be  the  demand  on  the  system.  FORM  uses  a  constrained  linear  optimization  technique  to  estimate  the 
probability  of  g  ^  0 . 
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The  main  benefit  of  using  probabilistic  reliability  techniques  such  as  FORM  to  predict  product 
performance  is  that  they  enable  the  engineer  to  determine  the  sensitivity  of  the  probability  of  failure  to 
individual  random  variables.  Specifically,  if  probabilistic  methods  are  used,  the  sensitivity  of  the  probability 
of  failure  to  changes  in  random  variable  mean  and  variance  can  be  found  in  closed  form  [19]. 

The  utilization  of  probabilistic  methods  to  estimate  product  performance  and  reliability  has  several 
benefits; 

1 .  Design  performance  sensitivities  to  changes  in  mean,  variance  and  truncation  of  variables  can  be 
estimated  using  closed-form  relationships. 

2.  Individual  and  system  level  estimates  of  performance  and  reliability  can  be  found  simultaneously. 

3.  Product  sensitivities  to  performance,  reliability  and  cost  can  be  linked  and  optimization  methods  can  be 
employed  to  determine  the  optimal  product  configuration. 


5.2. 1  Response  surface  modeling  of  the  limit  state 

The  limit  state  function  in  Eq.  (4.3)  is  not  available  as  a  closed-form 
function  of  the  basic  random  variables.  It  can  only  be  computed  through  die 
ANSYS  finite  element  stress  analysis;  thus  it  is  an  implicit  function  of  the 
random  variables.  In  such  a  case,  the  response  surface  approach,  as  shown 
in  Figure  10,  can  be  used  to  develop  an  approximate  closed-form 
expression  of  the  limit  state  function,  and  then  the  first-order  reliability 
method  (FORM)  can  be  used  to  estimate  the  delamination  probability.  The 
response  surface  approach  consists  of  two  steps: 

1 .  Design  of  experiments,  i.e.,  the  choice  of  different  sets  of  values  of  the 
random  variables  to  compute  the  corresponding  values  of  the  response  g, 
and 

2.  Regression  analysis,  to  construct  the  approximate  closed-form  function. 

An  initial  regression  analysis  (screening  analysis)  including  all  seven 
random  variables  revealed  that  only  three  out  of  the  seven  random 
variables,  namely,  the  bending  angle  {G),  magnitude  of  the  axial  load  (P), 
and  the  material  property  En,  had  any  significant  effect  on  the  limit  state 
function.  Therefore,  for  further  analysis  only  three  of  the  above  mentioned 
random  variables  were  included. 

A  Central  Composite  Design  (CCD)  (see  Section  6.3)  scheme  was 
adopted  to  fit  the  response  surface.  For  small  number  of  random  variables, 
such  a  design  scheme  is  most  efficient  to  explore  the  possibility  of  second- 
order  effects  and  first  order  interaction  terms  among  the  random  variables. 
A  total  of  fifteen  sets  of  experiments  (non-linear  finite  element  analyses) 
were  performed  to  develop  the  response  surface.  The  data  sets  for  each  experiment  were  gathered  by 
selecting  the  value  of  each  random  variable  at  tiu-ee  levels,  namely,  mean,  mean  +  2  standard  deviations  and, 
mean-  2  standard  deviations. 

Table  3  shows  the  array  of  15  experiments.  A  least  squares  regression  of  the  data  produced  the  equation: 
<Jtnax  =  1.006  +  8.402 X  10‘®^E:„  +0.138P  +  0.191^+5.432 x  10'^ 

+1.03 1x10-^^^+  2.359  X 10”^  PO 


Design  Matrix 
(input) 

X 


Figure  10:  Response 
surface  method. 
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As  expected,  the  regression  model  (Eq.(4.4))  showed  strong  effect  of  the  angle  and  magnitude  of  the  load. 
Among  the  material  property  variables,  Eji  was  also  found  to  have  some  significance.  Substituting  Eq  (4.4) 
into  Eq  (4.3)  and  maintaining  consistency  in  units  the  final  limit  state  equation  is: 


g  =  G„i,  - 175.344(1.006  +  8.402  x  10^’£:„  +0.138i’  +  0.1916»  + 5.432  x  10'"  P" 
+1.031  X 10'^  0^  +2.359  X  10'^  Pi?) 


(4.5) 


Table  3:  Response  surface  analysis  data 


No. 

P„(Msi) 

P(Lb) 

^Deg) 

Gmax(in-lb/in'‘) 

1 

6.90 

30800 

12 

1.014553 

2 

7.08 

30800 

12 

0.998717 

3 

6.72 

30800 

12 

1.031467 

4 

6.90 

36960 

12 

1.273955 

5 

6.90 

24640 

12 

0.668086 

6 

6.90 

30800 

15.333 

1.434791 

7 

6.90 

30800 

8.667 

0.633156 

8 

7.08 

36960 

15.333 

9 

7.08 

36960 

8.667 

0.817964 

po^ 

7.08 

24640 

15.333 

0.99471 

11 

7.08  ^ 

24640 

8.667 

0.468271 

12 

6.72 

36960 

15.333 

1.774213 

13 

6.72 

36960 

8.667 

0.832335 

14 

6.72 

24640 

15.333 

1.09308 

15 

6.72 

24640 

8.667 

0.483494 

5.2.2  Reliability  analysis 

The  objective  of  the  Phase  I  research  is  to  predict  the  probability  of  delamination  onset  in  a  helicopter 
rotor  hub  at  a  particular  number  of  loading  cycles.  The  steps  involved  in  estimating  the  probability  are: 

1 .  Choose  the  number  of  loading  cycles  (N). 

2.  From  Eq.  (4.2)  compute  the  mean  value  of  Gcrii  associated  with  N. 

3.  With  the  limit  state  equation  as  defined  in  Eq.  (4.5) ,  perform  a  FORM  analysis  to  estimate  the 
probability  of  delamination  onset  and  the  sensitivity  of  the  random  variables  to  the  estimated  probability. 

Figure  1 1  shows  the  estimated  probability  of  delamination  onset  at  different  values  of  the  loading  cycles. 
The  mean  predicted  life  is  approximately  15000  cycles.  Figure  12  shows  the  CDF  plotted  on  lognormal 
paper.  The  linear  form  of  the  CDF  in  Figure  12  indicates  that  the  fatigue  life  distribution  is  lognormal.  The 
sensitivity  results  of  the  FORM  analysis  are  shown  in  Table  4.  It  can  be  seen  from  the  results  that  the  bending 
angle  6,  the  axial  load  P,  and  Gcru  are  the  most  sensitive  variables. 
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Table  4  :Sensitivity  results 


Logl0(A) 

Gcrit 

P 

0 

1 

-0.532 

0.457 

0.712 

2 

-0.564 

0.452 

0.690 

3 

-0.599 

0.447 

4 

-0.0256 

0.442 

0.631 

5 

-0.0272 

0.437 

0.592 

6 

-0.713 

-0.0287 

0.436 

0.548 

5.3  Focus  Area  3:  Additional  Reliability  Algorithms 

The  composite  design  framework  developed  in  the  Phase  I  effort  allow  for  two  basic  types  of  reliability 
analysis;  response  surface  methods  and  mean  value  methods.  The  response  surface  method  is  best  suited  for 
the  construction  of  the  complete  CDF.  The  complete  CDF  analysis  is  useful  when  the  reliability  is  desired 
over  a  wide  range  of  cyclic  lives.  For  example,  the  designer  may  be  interested  in  knowing  the  reliability  for 
lives  ranging  from  1000  to  100,000  cycles.  In  this  case,  multiple  finite  element  analyses  are  performed  to 
create  a  global  response  surface  that  is  approximate  throughout  the  cyclic  life  range  of  interest.  The  response 
surface  module  includes  a  design  of  experiments  routine  to  determine  the  matrix  of  input  parameter  values 
based  on  the  desired  fidelity  of  the  response  surface.  The  response  surface  method  was  used  in  the 
demonstration  problem  to  determine  the  reliability  of  the  tapered  composite  as  discussed  in  Section  5.  This 
section  discusses  the  development  of  mean  value-based  reliability  algorithms  that  have  been  added  to  the 
overall  analysis  framework. 

Mean  value  methods  restricts  the  finite  element  analyses  to  the  region  of  the  most  probable  failure  point. 
This  technique  is  useful  when  the  reliability  at  a  single  value  of  life  is  of  interest.  For  example,  the  designer 
may  be  interested  in  knowing  the  reliability  only  at  a  design  requirement  of  10,000  cycles.  In  this  case,  it  is 
more  efficient  to  perform  multiple  finite  element  analyses  that  are  used  to  calculate  the  local  response  in  the 
vicinity  of  10,000  cycles.  Two  different  routines  were  developed  to  perform  the  local  reliability  analysis;  one 
based  on  Advanced  Mean  Value  methods  (AMV)  and  the  other  based  on  direct  first  order  reliability 
techniques. 

J.  3. 1  Mem  value  methods 

Mean  value  methods  are  particularly  suitable  when  a  closed-form  expression  of  the  limit  state  is  not 
directly  available,  as  is  the  case  when  finite  element  analysis  is  used  to  model  the  structural  response.  The 
reliability  computation  algorithm  developed  for  the  Phase  I  effort  constructs  a  linear  approximation  to  the 
limit  state  in  the  equivalent  standard  normal  space.  Mean  value  methods  works  differently,  as  follows. 

1 .  A  first-order  Taylor  series  approximation  is  constructed  at  the  mean  values  of  the  random  variables — in 
the  original  space— based  on  perturbation-based  sensitivity  analysis  using  the  finite  element  analysis  or 
other  implicit  response  models.  The  two-parameter  scheme  is  used  to  transform  this  closed-form 
approximation  and  the  random  variables  to  the  equivalent  standard  normal  space.  The  Rackwitz-Fiessler 
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iteration  formula  is  used  to  estimate  the  MPP,  and  a  first-order  estimate  of  the  failure  probability  is 
obtained  similar  to  FORM.  This  is  referred  to  as  the  mean  value  first-order  (MVFO)  estimate. 

2.  The  MVFO  estimate  is  refined  using  advanced  mean  value  (AMV)  analysis.  This  is  simply  a 
deterministic  analysis  of  the  system  at  the  MPP,  to  evaluate  the  g-function.  If  the  MVFO  estimate  were 
accurate,  the  g-function  would  be  exactly  zero.  If  the  MVFO  estimate  were  not  accurate,  then  a  different 
value  would  be  obtained  for  the  g-function.  The  MPP  and  the  probability  estimate  identified  by  the 
MVFO  are  now  assumed  to  correspond  to  this  new  value  of  the  g-function. 

3.  The  above  two  steps  are  repeated  for  different  values  of  the  g-function.  This  results  in  the  construction  of 
the  CDF  of  the  g-llinction. 


Mean  value  methods  are  a  practical  alternative  to  FORM,  if  the  deterministic  analysis  of  the  system  is 
expensive.  By  simply  combining  the  information  from  the  MVFO  step  and  one  additional  deterministic 
analysis,  one  is  able  to  obtain  a  substantially  improved  estimate  of  the  failure  probability. 

There  are  optimization  algorithms,  such  as  BEGS,  which  combine  the  information  in  the  iteration  steps  to 
perform  an  accurate  search.  Such  algorithms  are  especially  useful  when  the  Rackwitz-Fiessler  algorithm  fails 
to  converge.  However,  the  programming  effort  and  the  memory  storage  requirements  in  these  methods,  plus 
the  marginal  improvement  in  accuracy  over  the  Rackwitz-Fiessler  method  have  hindered  their  widespread 
application  in  structural  reliability  studies.  The  AMV,  with  its  minimal  one  step  combination,  is  a  practical 
approach  in  this  direction. 

The  AMV  estimate  can  be  further  improved  with  more  iteration.  In  the  AMV  +  method,  sensitivity 
analysis  is  performed  after  the  AMV  step,  and  steps  1  and  2  of  the  above  procedure  are  repeated.  In  a  similar 
manner,  further  iterations  can  be  done  to  produce  AMV  +  +,  AMV  +  +  +,  etc.,  estimates.  In  each  case,  note 
that  only  the  first  iteration  (MVFO)  and  the  latest  deterministic  analysis  are  combined.  It  is  important  to  note 
here  that  the  sensitivity  information  in  the  AMV  is  computed  only  in  the  first  iteration,  i.e.,  at  the  MVFO 
step,  which  is  at  the  mean  values  of  the  random  variables.  If  the  sensitivity  values  are  needed  at  the  MPP, 
then  full  sensitivity  analysis  needs  to  be  performed. 

The  AMV  method  achieves  accuracy  through  a  quadratic  approximation  of  the  g  function  at  the  very 
beginning  (at  the  mean  values  of  the  random  variables).  That  is,  the  first-order  derivatives  can  be  calculated  at 
the  beginning  using  finite  difference  or  perturbation  analysis  to  have  a  second-order  Taylor  series 
approximation 


g(X)  «  g(x*)  +  2(x*)  .  (x,.  - X.)  +  (Xi  - x'f 

/-I  '  1=1  ^ 

-X*)  +  ^6,(x, -x'f 


i=l 


1=1 


(5.1) 


Most  reliability  codes  available  today  linearize  the  above  expression.  However,  the  use  of  a  quadratic 
expression  for  the  g  function  in  Eq.  (5.1)  make  the  MPP  and  sensitivity  estimates  after  the  MVFO  estimate 
more  accurate.  Notice,  however,  that  the  quadratic  approximation  of  Eq.  (5.1)  does  not  Involve  mixed  terms, 
and  the  effects  of  the  random  variables  on  the  g-function  are  considered  separately.  This  introduces  some 
approximation. 


5.3.2  MVFO/AMV software  program 


The  MVFO  and  AMV  methods  described  above  were  successfully  transformed  into  a  user-friendly 
com]iuter  program  This  section  gives  an  example  of  how  to  use  the  method  with  the  help  of  an  example  and 
screen  dumps  from  the  program.  The  mathematical  details  of  the  example  are  given  in  Section  6.2.3. 
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Example: 


Assume  we  have  a  cantilever  bean  of  length  L,  loaded  at  the  free  end  with  a  force  w.  We  wish  to  find  the 
vertical  displacement  at  the  free  end  of  the  beam.  We  know  that  the  load  w,  stiffness  E,  and  the  moment  of 
inertia/,  are  all  normally  distributed  random  variables,  with  mean  and  standard  deviation  as  follows: 

iV(5000,1000)  lbs 

I « iV(120,l0)  in 

£«#(30xl0®,0.3xl0®)psi 


/«iV(1200,150)in^ 

We  know  from  mechanics  that  the  vertical  displacement  d  of  the  free  end  of  a  beam  is 

wl} 


d  = 


3EI 


(5.2) 


Let  us  assume,  for  this  example,  we  could  not  determine  the  displacement  in  closed  form  and  had  to  use  a 
simple  finite  element  model  to  determine  the  end  displacement.  The  details  of  each  step  involved  in 
performing  a  Mean  Value  First  Order  (MVFO)  Analysis  for  the  above  example  are  enumerated  below. 

MVFO  Analysis 


Step  1.  The  first  step  is  to  input  all  of  the  above  information  into  the  computer  program.  Figure  13  and  Figure 
14  show  the  input  screens  for  this  step.  It  should  be  noted  that  random  variables  xl ,  x2,  x3  and  x4  refer  to 
variables  w,  /,  E  and  /  respectively. 

Step  2.  The  next  step  in  a  MVFO  analysis  is  performing  a  perturbation  analysis  about  the  mean  values  of  the 
input  random  variables  (w,  L,  E,  1)  using  the  finite  element  model.  The  perturbation  analysis  is  performed  to 
ck 

determine  the  gradient  —  of  the  response  to  each  of  the  random  variables.  In  this  example  z  represents 

ac, 

vertical  tip  displacement  d,  of  the  cantilever  beam.  The  perturbation  analysis  involved  changing  each  variable 
by  one-tenth  of  its  standard  deviation  about  its  mean  value  and  rerunning  the  finite  element  model  to  create 
the  array  of  five  finite  element  runs  as  shown  below.  The  results  of  the  perturbation  analysis  are  shown  in 
Figure  15. 

Step  3.  Based  on  the  perturbation  analysis  of  Step  2,  a  first  order  Taylor  Series  expansion  is  performed  about 
the  mean  values  of  the  random  variables  expressing  the  vertical  tip  displacement  cf,  as  a  function  of  the  input 
variables.  Figure  15  also  shows  this  Taylor  Series  expansion  expression. 

Step  4.  Having  defined  the  Taylor  Series  expansion,  the  next  step  is  to  define  the  limit  state  equation  for  first 
order  reliability  analysis.  The  limit  state  equation  is  defined  as: 

g  =  d^-dapprox  (5-3) 

where  is  the  "level  of  displacement"  for  which  the  reliability  analysis  is  performed  and  dt^x  is  the 
Taylor  Series  expansion.  In  other  words,  the  reliability  analysis  will  find  the  probability  that  the  beam 
displacement  greater  than  d^ .  In  the  program  d^  is  referred  to  as  "Z-Level"  (Figure  16).  The  user 

has  a  choice  (a)  either  suggesting  a  range  of  Z-Ievels  {d^  values)  based  on  a  constant  increment  value 
(Automatic  Z-Level  Input)  or  (b)  if  the  user  is  only  interested  in  specific  Z-levels,  he/she  can  input  those 
specific  values  individually  (Manual  Z-Level  Input). 
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Figure  13:  AMV  input  screen  -1. 


Figure  14:  AMV  input  screen  -  2. 
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Figure  15:  Perturbation  analysis 


Figure  16:  Z-level  input. 
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Based  on  a  range  of  values  of  do,  the  probability  that  g  <  0  is  determined.  Figure  1 7  lists  the  range  of 
values  of  do  considered  along  with  the  corresponding  reliability  indices  (P)  and  the  most  probable  values 
(MPP)  of  the  input  random  variables  (w,  L,  E,  I)  for  which  the  relationship  g  <  0  was  satisfied.  Figure  17  also 
plots  t/o  vs.  P  (z  vs.  beta)  for  the  MVFO  analysis,  which  represents  the  cumulative  distribution  function 
(CDF)  of  the  vertical  tip  displacement  d  of  the  beam.  Notice  that  the  MVFO  curve  is  linear.  This  is  because 
the  random  variables  are  normal  and  a  first-order  approximation  to  the  beam  displacement  response  (equation 
(5.3))  was  used  in  the  analysis.  However,  we  have  no  reason  to  believe  that  the  true  displacement  response  of 
the  beam  is  linear.  If  we  wish  to  develop  a  better  approximation  to  the  true  CDF,  we  perform  an  advanced 
mean  value  analysis. 

AMV  Analysis 


AMV  is  a  continuation  of  the  MVFO  analysis  in  which  we  use  the  MPP  from  the  MVFO  analysis  and 
the  finite  element  model  to  further  refine  the  reliability  analysis. 

Step  5.  The  vertical  tip  displacement  based  on  the  finite  element  analysis  using  the  MPP  values  of  the  random 
variables  is  evaluated.  These  values  are  listed  in  the  last  column  (z_org)  of  MVFO  Results  table  in  Figure  17. 
Notice  that  the  value  of  the  tip  displacement  from  the  finite  element  analysis  (zorg)  is  not  the  same  as 
(z-level).  The  value  of  (z-level)  is  equivalent  to  the  value  obtained  by  using  Taylor  Series  approximation 
at  the  MPP.  The  displacement  value  at  each  P  is  shifted  or  "moved"  from  Jg  to  dp^^  as  show  by  the 
nonlinear  AMV  curve  in  Figure  17.  The  AMV  curve  is  the  CDF  that  represents  a  refinement  to  the  MVFO 
CDF  using  just  one  additional  finite  element  run. 
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Figure  17:  MVFO  results. 
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5.3.3  Direct  FORM 


Direct  FORM  is  an  extension  of  the  MVFO  methods.  In  the  direct  FORM  method,  we  do  not  construct  a 
closed-form  approximation  to  the  limit  state.  The  finite  element  analysis  is  used  to  compute  the  value  and  the 
gradient  of  the  limit  state  at  each  iteration  of  FORM  is  used  to  search  for  the  most  probable  point  (MPP).  The 
number  of  iterations  is  based  on  the  convergence  of  and  the  MPP.  Thus,  AMV  can  be  considered  a  subset 
of  direct  FORM. 

To  implement  the  direct  FORM  method,  a  C  program  has  been  written  to  directly  integrate  finite  element 
analysis  and  reliability  analysis.  The  program  does  the  following  steps: 

1 .  Read  the  mean  value  and  standard  deviation  of  each  random  variable  as  input  data. 

2.  Call  the  laminate  code  to  calculate  the  material  properties  of  laminate  model. 

3 .  Generate  the  input  file  for  calling  the  finite  element  software  as  a  subroutine  to  perform  the  structural 
analysis  of  the  rotor  hub  component. 

4.  Read  the  finite  element  analysis  output  file  and  calculate  the  limit  state  function  g(x)  and  its 
derivatives  with  respect  to  the  random  variables. 

5.  Use  FORM  to  compute  the  next  iteration  point  in  the  search  for  the  MPP. 

Steps  2  to  5  are  repeated  at  each  iteration  until  the  program  converges  to  the  MPP.  This  provides  the  estimate 
of  the  reliability  index  p  and  the  delamination  probability  of  rotor  hub  component. 

The  direct  FORM  program  is  re-run  for  different  number  of  cycles  of  fatigue  life,  and  the  probability  of 
failure  vs.  number  of  cycles  is  computed.  For  this  problem,  the  results  of  direct  FORM  are  observed  to  be 
within  10%  of  the  results  with  the  response  surface  approach. 

6  Theoretical  Background  for  Reliability-Based  Design 
6.1  Interference  Theory 

Reliability-based  design  stems  from  the  concepts  of  interference  theory.  Interference  theory  is  used  to 
find  the  probability  that  the  resistance  (r)  of  a  component  or  system  is  less  that  the  applied  load  (/).  The 
simple  mathematical  function  that  describes  the  performance  of  the  design  is 

g  =  r-l  (6.1) 

which  is  called  the  response  or  performance  function.  The  regime  in  which  g  =  r-l>0  is  the  safe  state,  the 
regime  in  which  g  =  r  -  /  <  0  is  the  failure  state,  and  the  state  g  =  r-l  =  0  separates  the  safe  state  from  the 
failure  state  and  is  called  the  limit-state. 

The  term  interference  theory  comes  from  the  graphical  representation  of  the  PDF  for  the  resistance  and  load 
as  shown  show  in  Figure  18.  The  failure  state  can  be  associated  with  the  shaded  region.  This  region  is 
associated  with  the  probability  of  failure  or  the  probability  that  the  resistance  is  less  than  the  load, 

P(r<l)  =  P{r-l<0)  =  P(g<0) 

which  is  the  probability  that  the  resistance  random  variable  f  is  less  than  a  realization  of  the  load  /  multiplied 

A, 

by  the  probability  that  the  load  random  variable  /  is  equal  to  the  realization  I  for  all  realizations  of  the  load, 

P(r<l)  =  Y,P{r<l)Pil<l)  (6.2) 

all/ 

We  know  that  the  probability  that  r  is  less  than  /  is  equal  to  the  CDF  of  r  at  L 
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Figure  18:  Resistance  load  interference  diagram. 

P{f<l)  =  FXl)  ^  (6.3) 

and  the  probability  that  /  is  equal  to  /  is  the  probability  that  /  is  between  /  and  a  small  increment  A/ . 

P(/ </)  =  ?(/</ </  +  A/) 

P(i<l)  =  Ml)M 

Substituting  (6.3)  and  (6.4)  into  (6.2)  yields 

Pif<i)  =  '£K(i)MiW  (6.5) 

all/ 

The  expressions  in  (6.5)  are  represented  graphically  in  Figure  19. 

As  A/  becomes  infinitesimally  small  and  approaches  dl,  the  probability  of  failure  can  be  expressed  as 


(6.4) 


/>(r  </)  =  />,(/)/, (/)<//  (6.6) 

Alternatively,  we  can  consider  the  probability  of  failure  as  the  probability  that  the  load  is  greater  than  the 
resistance 


P(,t>e)  =  ^[\-F,(t)]fXr)dr  (6.7) 

For  particular  statistical  distributions  of  resistance  and  load,  direct  solutions  to  (6.6)  (or  (6.7))  are  available 
[Kapur  and  Lamberson,  1977].  If  f  and  I  are  independent  normal  (Gaussian)  random  variables,  then 
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(6.8) 


r~N{n„a,) 

I  ~N{n„a,) 
rFrmWl  =  ^{-P)  (6.9) 

wO 

p=  (6.10) 

Vo-;  +  cr/ 

where  O  is  the  standard  normal  CDF  and  P  is  referred  to  as  the  reliability  index  or  safety  index. 

Recall  that  for  linear  functions  of  independent  random  variables,  the  mean  value  of  the  function  is  equal 
to  the  function  evaluated  at  the  mean  values  of  the  individual  variables  and  therefore  the  mean  value  of 
g  =  r-l  is 


Mg=Fr-  Ml 

and  the  variance  of  g  is  the  square  root  of  the  sum  of  the  squares  of  the  variance  of  the  individual  random 
variables 

V2  2 

(J,+£7, 

The  safety  index  (6.10)  can  be  written  as 


M.-Mi  Mg 


+  0-/ 


(6.11) 


which  is  the  inverse  of  the  COV  of  g. 

If  r  and  /  are  independent  lognormal  random  variables,  then 

r  ~  LN{E{f),StdDev{r)) 

I  ~  LN{E{l),StdDev{l)) 
'StdDev{xf 


0-1  =  In 


E{x) 


-  +  1 


Mx=^^E{x)-^(tI 

\yri.i)mdi=^(.-P) 


P  = 


Mr- Ml 


(6.12) 


where  E(x)  is  the  expected  value  and  StdDev{x)  is  the  standard  deviation.  If  r  and  /  are  independent 
exponential  random  variables,  then 

E{r)  =  ^ 

1 


£(/)  = 


A, 

A,  E(f) 


^r-^l  Eif)  +  E(l) 
If  r  and  /  are  independent  and  r  is  normal  and  I  is  exponential,  then 

r  ~ 

i~EXPO{X,) 


(6.13) 
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Figure  20:  Interference  diagram  showing  increase  in  interference  area  with  decrease  in 

safety  measures. 


\  ^rj 


-exp 


i-d 


z"  „  32^2  ^ 

Af  r 


O'r 


(6.14) 


If  r  and  /  are  independent  and  r  is  exponential  and  /  is  normal,  then 

f~EXPOiA,) 


p;(/)/,(/)rf/=l-<D 


V  o-zy 


+  exp 


1 

1 

T 

e 

f  ^2  2 

JL 

^  ^1  J. 

(6.15) 


The  probability  of  failure  depends  on  the  mean  (or  average)  value  of  the  resistance  and  load.  As  the 
difference  between  the  mean  resistance  and  mean  load  increases,  the  probability  of  failure  decreases.  If  the 
difference  between  the  mean  resistance  and  mean  load  decreases,  the  probability  of  failure  increases.  This  is 
shown  graphically  in  Figure  20.  Note  that  the  ratio  of  the  mean  resistance  to  the  mean  load  is  equivalent  to 
the  safety  factor  used  in  deterministic  analysis. 


safety  factor  = 


(6.16) 


The  difference  between  the  mean  resistance  and  the  mean  load  is  the  safety  margin  used  in  deterministic 
analysis. 


safety  margin  =  /j,  (6.17) 


As  the  safety  factor  or  safety  margin  decrease,  the  probability  of  failure  increases. 

The  probability  of  failure  also  depends  on  the  scatter  or  degree  or  dispersion  of  the  resistance  and  load.  If 
either  the  scatter  in  resistance  or  the  scatter  in  load  increases,  the  probability  of  failure  increases.  Figure  2 1 
shows  that  as  the  scatter  increases,  the  size  of  the  interference  region  increase.  The  change  in  probability  of 
failure  with  changes  in  the  mean  values  and  scatter  of  the  random  variables  can  also  be  observed  by 
considering  the  magnitude  of  the  safety  index  /?  as  shown  in  (6.10). 

We  have  shown  two  methods  for  decreasing  the  probability  of  failure  of  a  design.  One  method  is  to 
increase  the  relative  distance  between  the  mean  resistance  and  the  mean  load.  This  is  the  result  of  traditional 
design  improvements  which  is  equivalent  to  increasing  the  safety  factor.  The  other  method,  which  traditional 
design  improvement  does  not  address,  is  to  decrease  the  scatter  in  the  resistance  and/or  load  variables. 
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Figure  21:  Interference  diagram  showing  increase  in  interference  area  with  increase  in  scatter. 


Addressing  the  scatter  in  the  design  parameters  gives  the  design  a  whole  different  method  to  improve  the 
reliability,  efficiency,  and  competitiveness  of  a  design.  Increasing  the  safety  factor  has  been  the  most  popular 
design  improvement  method  because  methods  addressing  scatter  have  not  been  widely  available  except  for 
the  simple  r-l  performance  function  and  a  limited  variety  of  statistical  distributions.  Tabulated  numerical 
solutions  to  (6.6)  for  other  distribution  type  (most  notably  Weibull)  are  available  (see  Kapur  and  Lamberson, 
1977).  However,  these  solutions  are  still  based  on  the  simple  r-l  performance  function. 

Let  us  consider  a  simple  yet  realistic  design  of  a  cantilever  beam.  Assume  the  beam  is  subjected  to  a 
point  load  w  at  the  free  end.  The  design  criteria  specifies  that  the  displacement  d  at  the  free  end  does  not 
exceed  the  a  specific  value  denoted  by  z.  From  mechanics  of  materials,  we  know  a  performance  function 
describing  the  end  displacement  is 


d  = 


wL^ 

3EI 


(6.18) 


where  w  is  the  load,  L  is  the  length  of  the  beam,  E  is  Young’s  modulus  of  the  beam  material,  and  I  is  the 
cross  sectional  moment  of  the  beam.  The  limit  state  that  separates  the  safe  region  from  the  failure  region  is 

g  =  z-d  =  z-  —  =  0  (6.19) 

3EI 


Assume  the  deflection  criteria  z  is  deterministic  and  the  load,  length,  modulus,  and  moment  are  all  random 
variables  with  the  following  distributions; 


L~LN{A„C,) 

I~LN{A.„C,) 


E~WB{b„tjE,r,) 


The  probability  of  failure  is  formulated  as 


Pf=P\ 


z<- 


3EI 


(6.20) 


We  rewrite  (6.20)  so  that  we  have  random  variables  on  both  sides  of  the  inequality 
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(6.21) 


PT/'  i-^^P  (lognormal)  (6.24) 

=  ~T  —y^'  (lognormal)  (6.25) 

\27r^iJ  _  ■^  \  ^1  J 

h  f  F  -  v^”'  r  ( F  -  ' 

f^{E)  =  -^  - ^  exp - ^  (Weibull)  (6.26) 

IeK  rjE  )  [  I  »?£  J 

Although  the  performance  function  describing  the  displacement  of  the  beam  (6.1 8)  is  simple,  substituting 
(6.23)  through  (6.26)  into  (6.22)  creates  a  multiple  integral  that  is  difficult,  if  not  impossible  to  solve.  Using 
direct  integration  to  find  the  probability  of  failure  for  most  engineering  designs  is  not  practical. 

6.2  Second  Moment  Formulation 

Due  to  the  impracticality  of  solving  the  multiple  integral  associated  with  the  interference  region,  second 
moment  formulations  have  been  developed  to  approximate  the  probability  of  failure.  The  term  second 
moment  formulation  refers  to  the  transformation  of  the  random  variable  statistics  into  equivalent  standard 
normal  variables  described  only  by  the  first  and  second  moments,  p  and  irrespectively.  We  will  develop  the 
second  moment  formulation  based  on  the  simple  r-l  performance  function  with  both  f  and  /  as 
independent  Gaussian  random  variables.  The  formulation  will  then  be  generalized  for  more  complicated 
performance  functions  and  various  types  of  statistical  distributions. 

6. 2. 1  Linear  response  functions 

If  g  =  r  -  /  then  we  observe  that 

(g  >  0)  -)•  (r  >  /)  is  the  safe  state 
(g  <  0)  (r  <  /)  is  the  failure  state 
(g  =  0)  ->  (r  =  /)  is  the  limit  state 

If  we  transform  r  ->  A(^,.,cr,.)  and  /  ->  N{p„(T,)  into  standard  normal  variables  f'  ->  A(1,0)  and 
/ '  ^  iV(l,0) 

f'  =  f  =  o-,f '  +  //,  (6.27) 

(T^ 
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Figure  22:  Linear  limit  state  equation  in  standard  normal  space. 
7'  =  ^^:^^/ =<T, /'  +  //,  (6.28) 

o-/ 

The  limit  state  equation  can  now  be  written  as 


g=:0  =  f-/ +  /'-//,  (6.29) 

If  we  plot  (6.29)  in  standard  normal  space  ( considering  the  axis),  we  get  a  straight  line  as  shown  in 

Figure  22. 

We  know  that  the  minimum  (perpendicular)  distance  S  of  any  line 

ax  +  Z)y  +  c  =  0  (6.30) 

to  the  origin  is 


(6.31) 


If  we  put  (6.29)  in  the  form  of  (6.30) 

g  =  [a,]?'  +  {-a,)V  +  (//,-//,)  =  0 
then  the  minimum  distance  S  to  the  origin  is 


(6.32) 


5  = 


■  Mr- Ml 


+  a] 


(6.33) 


Comparing  (6.33)  to  (6.10)  we  notice  that  5  is  equal  to  the  safety  index  P  and  thus  the  probability  of  failure 
for  the  performance  function  is 

Pj  =  cD(-<J)  (6.34) 


The  point  on  the  limit  state  line  which  is  closest  to  the  origin  is  called  the  most  probable  point  (MPP) 
designated  j .  Of  all  of  the  values  of  which  cause  failure,  the  MPP  value  can  be  thought  of,  in 

some  approximate  sense,  as  the  value  of  which  are  most  likely  to  occur.  The  MPP  is  an  important 

concept  in  the  second  moment  formulation  because  the  distance  from  the  MPP  to  the  origin  is  the  safety 
index.  Also,  approximation  methods  that  allow  us  to  consider  complex  performance  functions  must  make  use 
of  the  MPP. 
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The  x'j  coordinate  of  the  MPP  can  be  found  using  analytical  geometry 

x';=-a]p  (6.35) 

where  or*  is  the  cosine  of  the  angle  between  the  P  vector  and  the  XI  axis 

^1 


or,  = 


.=1  K^'iJ 


(6.36) 


(We  will  drive  (6.36)  later.)  If  we  apply  the  above  formulations  to  the  simple  linear  performance  function 

A 

g  =  r-l  with  both  r  and  /  as  independent  Gaussian  random  variables,  the  limit  state  function  transformed 
into  standard  normal  variables  is 

g  =  {cr,)r'  +  {-<J,)V  +  (//,-;/,)  =  0  (6.37) 

with  safety  index 

The  MPP  in  standard  normal  space  is  designated  (/■'*,/'*)  where 

(/■'*,/'*)  =  (-a 

dr' 


cP 


rV 

')  U'J 


a,  =■ 


a' 


a^'J  [a' 


Differentiating  (6.37)  with  respect  to  f'  and  /'  respectively 


Therefore, 


a' 


=  cr. 


ar  = 


=  0, 


a.  = 


-cr, 


and  the  MPP  in  standard  normal  space  is 


-GrP  a,P 


Vo-r+O-/  Vo’r+O'? 
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The  direction  cosines  a*  are  also  called  the  variable  importance  factors  (VIF).  They  provide  a  relative 
measure  of  the  sensitivity  of  the  reliability  to  the  random  variable  Xj .  The  sign  of  a]  indicates  the  influence 
of  the  variable  on  the  reliability.  The  positive  sign  on  a*  indicates  that  as  the  resistance  increases,  the 
reliability  increases.  The  negative  sign  on  a]  indicates  that  as  the  load  increases,  the  reliability  decreases. 

The  MPP  in  real  space  is 


(r-/) 

V  J 

Numerical  Example  6.1:  Assume  we  are  interested  in  analyzing  an  automatic  assembly  process  in 
which  3  steel  sheets  are  stacked  on  top  or  each  other  and  inserted  into  a  plastic  clip  as  shown  in 
Figure  23. 


-o.B  (jiB  '' 

/  2  2~Br  yr . f-Bl 

<T,  <7/ 


Plastic  Clip 


Steel  Sheets 


Figure  23:  Steel  sheets  insert  into  the  plastic  clip. 


The  assembly  process  fails  when  the  stack  of  steel  sheets  is  too  thick  to  fit  into  the  plastic  clip.  The 
performance  function  is  defined  as 

g  =  x^-Zx,  (6.38) 

We  define  the  probability  of  failure  as  the  probability  that  the  clip  opening  is  less  than  the  stack 
thickness. 

Pj=P(xp<2x,)  (6.39) 

The  limit  state  equation  associated  with  the  (6.38)  is 

g  =  Q  =  Xp-2x,  (6.40) 

We  know  that  measurements  made  during  statistical  process  control  indicate  that  the  sheet 
thickness  and  the  clip  opening  have  the  following  statistical  distributions: 

Xs~N{p„<j,)  =  N{\  in,  0.1  in) 

Xp  ~  N{pp,ap)  =  N{A  in,  0.8  in) 

The  steel  sheets  and  the  plastic  clips  are  manufactured  by  completely  different  methods.  So.  it  is  fair 
to  assume  that  the  sheet  thickness  and  the  clip  opening  are  independent  random  variables. 
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To  find  the  probability  of  failure  we  transform  the  random  variables  into  standard  normal 

variables 


x'  =  — 


■Ms 


'Xs=<ysXs+Ms 


Substituting  (6.41)  into  (6.38)  yields 

g = 0 = (cTpji; + (-3<rjje; + -  3//,) 


(6.41) 


(6.42) 


Comparing  (6.42)  to  (6.32)  allow  us  to  solve  for  the  safety  index  /? ,  which  is  the  minimum  distance 
from  the  line  g  =  0  to  the  origin  as  in  (6.33). 


The  probability  of  failure  is 


Mp-^Ms  4-3 

"Vo.8'-0.3= 


1.17 


=(l)(-y3)  =  (I)(-1.17)  =  0.121 


Which  indicates  that  about  12%  of  the  assemblies  will  fail. 


It  is  interesting  to  note  that  the  deterministic  analysis  tell  us  that  we  have  a  safety  margin  of  1  inch 
and  a  safety  factor  of 

safety  factor  =  =  —  =  1.33 

^Ms  3 

If  we  wish  to  change  the  design  parameters  such  that  we  decrease  the  probability  of  failure,  it  is 
useful  to  calculate  the  VI  Fs  from  (6.36) 


0.8 


Vo.8^  +  0.3^ 


=  0.936 


a,  = 


-30-. 


-0.3 


+{2a,f  Vo.8^  +0.3^ 


=  -0.351 


The  VI  Fs  indicate  that  the  probability  of  failure  is  most  sensitive  to  changes  in  Xp .  Let’s  check  this 
out. 


Assume  that  through  an  improved  manufacturing  process,  the  plastic  clip  can  be  manufactured  with 
tighter  tolerance  such  that  the  standard  deviation  is  cut  in  half  and  the  steel  sheets  remain 
unchanged 

Xp~N{4,0A) 

Xs-Nim) 

The  new  probability  of  failure  is  found  as  follows: 

/3  =  —^z2=  =  2 
■s/o.4^  +0.3^ 

Pj  =0(-2)  =  0.0228 

Which  indicates  that  about  2%  of  the  assemblies  will  fail  which  is  a  substantial  improvement  over  the 
12%  probability  of  failure  for  the  original  assembly  procedure. 
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Now  let  us  assume  that  through  an  improved  manufacturing  process,  the  steel  sheets  can  be 
manufactured  with  tighter  tolerance  such  that  the  standard  deviation  is  cut  in  half  and  the  plastic  clip 
remains  unchanged 

x,~m0.8) 
a:,  ~  iV(l,0.05) 

The  new  probability  of  failure  is  found  as  follows: 

4-3 

=1.23 

voF+oTi^ 

Pj  =(I)(-1.23)  =  0.109 

Which  indicates  that  about  11%  of  the  assemblies  will  fail  which  is  not  a  substantial  improvement 
over  the  12%  probability  of  failure  for  the  original  assembly  procedure. 

Therefore,  it  would  be  wise  to  address  the  manufacturing  tolerances  of  the  plastic  clip  if  a  decrease 
in  the  number  of  assembly  failures  is  desired.  Deterministic  analysis  methods  do  not  allow  us  to 
assess  the  improvements  gained  through  the  tighter  manufacturing  tolerances.  Deterministic  design 
analysis  would  only  allow  us  to  change  the  mean  values  of  the  random  variables. 


We  can  extend  the  above  example  for  the  generalized  linear  performance  function  of  the  form 

g(x)  =  Oo  +  ^a^x,  (6.43) 

1  =  1 

The  corresponding  limit  state  equation  is 

n 


/=! 


If  the  random  variables  are  independent  Gaussian  variables  with  distribution  ,  we  can  transform 

them  into  standard  normal  variables 


O’, 


and  express  the  limit  state  equation  as 

n 

«o  +  +  Pi)  =  0  (6.45) 

1=1 

We  know  that  the  distance  S  from  any  surface 

Co+ZQt,=0  (6.46) 

to  the  origin  is 


(=1 


Cn 


Writing  (6.45)  in  the  form  of  (6.46)  yields 


(6.47) 


(6.48) 


1=1  y  ;=1 


The  safety  index  P  can  be  expressed  as  the  minimum  distance  S  by  comparing  (6.46),  (6.47)  and  (6.48)  as 
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6.2.2  General  performance  function 


So  far  we  have  only  considered  the  case  of  the  linear  performance  function.  Now  consider  the  case  of 
generalized  performance  function  of  the  form 

g(^)  =  g(A:|,jC2---,x„)  (6.50) 

where  x  is  the  vector  or  n  independent  Gaussian  random  variables.  If  we  replace  each  random  variable  x, 
by  its  equivalent  standard  normal  variable 

'  (7,  (6.51) 

x,=c7,x;  +  pi 

the  limit-state  equation  can  be  written  as 

g(Xi,X2---,X„)^g{cr^xl  +Pi,(72X2  +*“»)  =  0  (6.52) 

The  distance  D  from  any  point  x'  =  {x{,X2,---,x'„)  in  standard  normal  space  to  the  origin  of  standard  normal 
space  is 

D  =  -Jxl^  +X2^+---+xl,^  =  (6.53) 

The  MPP  of  the  general  performance  function  (which  is  the  point  on  the  limit-state  equation  (6.52)  closest  to 
the  origin)  can  be  found  using  optimization  techniques  by  minimizing  D  subject  to  the  constrain  that 
g(x)  =  0 . 

Lagrangian  multipliers  may  be  used  for  this  optimization  as 

L  =  D  +  2gix)  (6.54) 

or  in  scalar  form 

+  (5  55) 

+^s(x„) 

where  (6.55)  represents  a  set  of  «  equations  and  X/  and  x/  are  related  according  to  (6.51).  To  find  the 
minimum  D  we  minimize  L  as 

^  =  0  i  =  l,2,  -,tt  (6.56) 

Scl 

and 

f  =  0,6.57) 

which  yields  n  + 1  equations  with  n  + 1  unknowns.  The  solution  will  provide  x'  which  is  the  MPP  in 
standpd  normal  space  Jc'*  =  (x'*i  ,x’*2  . 

Substituting  (6.54)  into  (6.56)  yields 

-^  =  ^  +  2-^  =  0  1  =  1, 2, (6.58) 

axl  D  aci 

from  which 
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Figure  24:  Nonlinear  limit  state  equation  in  standard  normal  space. 


The  minimum  distance  p  is  found  when  D  is  evaluated  at  the  MPP. 


(6.65) 


where  the  subscript  *  indicates  that  the  expression  is  numerically  evaluated  at  the  MPP.  Later,  we  will  present 
numerical  examples  to  find  p  using  the  equations  outlined  above.  Before  we  do  that,  however,  we  should 
discuss  the  significance  of  P  for  the  nonlinear  performance  function. 

Remember  that  p  is  the  distance  from  the  origin  to  the  MPP  in  standard  normal  space  for  linear  limit 
state  equations  as  shown  previously  in  Figure  22  and  nonlinear  limit  state  equations  as  shown  in  Figure  24. 
Recall  that  in  section  6.2.1,  we  showed  that  p  is  the  safety  index  and  can  be  used  to  find  the  probability  of 
failure  for  a  linear  performance  function  as 

Pj=<^{-P)  (6.66) 


This  is  an  exact  equation  with  no  approximations.  However,  p  cannot  be  used  to  find  the  exact  probability  of 
failure  for  a  nonlinear  performance  function.  We  can  use  P  to  approximate  the  safety  index  if  we  consider 
the  following. 

Figure  25  compares  a  linear  limit  state  equation  to  a  nonlinear  limit  state  equation  with  the  same  MPP 
and  p .  A  line  along  p  forms  a  vector  that  is  perpendicular  to  the  line  form  by  the  linear  limit  state  equation 
gi{X)  =  0 .  Also,  the  line  formed  by  the  linear  limit  state  equation  is  tangent  to  the  line  formed  by  the 


Figure  25:  Linear  and  nonlinear  limit  state  equations  in  standard  normal  space. 
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nonlinear  limit  state  equation  gi^(X)  =  0  . 

We  know  from  (6.66)  that  /?  can  be  used  to  find  the  probability  of  failure  of  the  linear  limit  state.  The 
linear  limit  state  is  equivalent  to  a  linear  approximation  to  the  nonlinear  limit  state  at  the  MPP.  Therefore, 
can  be  used  in  (6.66)  to  find  the  probability  of  failure  for  the  first  order  approximation  of  the  high  order  limit 
state.  Using  f  in  such  a  way  is  the  basis  for  what  is  know  as  first  order  reliability  methods  or  FORM.  We 
shall  show  later  that  FORM  provides  good  approximations  for  the  probability  of  failure  for  most  nonlinear 
limit  states.  This  is  primarily  because  the  transformation  of  the  random  variables  to  standard  normal  variables 
often  serves  to  rectify  or  straighten  the  curve  formed  by  the  non-linear  limit  state  equation. 

An  alternate  and  much  simpler  way  of  deriving  the  results  presented  in  (6.65)  through  (6.76)  is  to  first 
perform  a  linear  approximation  to  the  general  limit  state  equation  and  then  derive  an  expression  for  f  [20]. 
We  can  expand  the  general  nonlinear  performance  function  at  the  MPP  in  a  Taylor  series  as 

g(x)  =  g(x’)  +  X(x,-x;)f^]  +HOT  (6.67) 

where  HOT  represents  the  higher  ordered  terms  in  the  series  expansion.  We  know  that  g^Jc* )  =  0  because  the 

MPP  lies  on  the  failure  surface  and  if  we  ignore  the  HOT,  then  we  can  write  the  first  ordered  approximate  to 
the  performance  function  as 

J* 

Recall  that  jc,  =  +  ju^  and  therefore  we  can  write 

X, -x' +  +Mi)  =  cr^(x'  +  x';)  (6.69) 

and 


^  3cl  1  (  ^ 


3c j  3c\  3c j  cr,  {^3cl 


(6.70) 


Substituting  (6.69)  and  (6.70)  into  (6.68)  yields 


f;  (6-71) 

1=1  \^i  j* 

Equation  (6.71)  is  represented  by  the  limit  state  equation  gfX)  =  0  in  Figure  25,  which  is  the  linear 
approximation  to  the  general  nonlinear  limit  state  equation  gj^iX)  -  0  also  shown  in  Figure  25.  Recall  that 
the  mean  value  of  the  standard  normal  variable  is  0,  thus,  the  mean  value  of  g(x)  is  approximated  by 


u  =-yjc'*  ^ 
;=1 


(6.72) 


Recall  that  the  variance  of  the  standard  normal  variable  is  1,  thus,  the  variance  of  g(x)  is  approximated  by 


n  f  A^\  f  A,  Y 


We  know  from  (6. 1 1 )  that  f  is  the  ratio  of  the  mean  value  to  the  standard  deviation  of  g(x) , 


which  is  the  same  as  (6.65). 
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Remembering  that  each  coordinate  of  the  MPP  can  be  expressed  as 

xf=-a]p  (6.74) 

and  replacing  Z)  with  P  in  (6.61)  yields, 

f  XU  ^ 

ly"  ^ 

V'  '  y 

Substituting  (6.74)  into  (6.5 1)  yields 

X*  =fJi-aiGiP 

Equations  (6.75)  and  (6.76)  can  be  used  in  a  numerical  algorithm  to  solve  for  the  MPP  as  follows  [20,21] 

1 .  Assume  coordinates  for  the  MPP  (usually  the  mean  values  of  the  random  variables  is  a  good  starting 
point). 

2.  Derive  each  of  the  derivatives  — . 

Sx'j 

3.  Evaluate  each  of  the  derivatives  at  the  assumed  coordinates  of  the  MPP. 

4.  Evaluate  each  of  the  direction  cosines  a]  according  to  (6.75). 

5.  Write  each  of  the  coordinates  of  the  MPP  x*  in  terms  of  p  according  to  (6.76). 

6.  Substitute  the  x*  ’s  found  above  into  the  limit  state  equation  g(x)  =  0  such  that  the  limit  state  equation  is 
in  terms  of  P . 

7.  Solve  for  p . 

8.  Find  the  updated  assumed  coordinates  of  the  MPP  x]  by  reevaluating  (6.76)  with  the  p  found  above. 

9.  Repeat  steps  3  through  8  until  convergence  is  obtained  on  p . 

There  are  many  other  FORM  algorithms  to  solve  for  /?  [22].  There  are  also  methods  available  to  find  the 
MPP  in  real  space  i.e.,  no  transformation  to  standard  normal  space  is  performed  [23].  There  are  also  many 
other  mathematical  algorithms  available  in  reliability  analysis  such  as  second  order  reliability  methods  or 
SORM  and  direct  numerical  solutions  to  the  integral  shown  in  (6.6).  But,  for  illustrative  purposes,  we  will 
use  the  simple  algorithm  outlined  above. 


(6.75) 


(6.76) 


Numerical  Example  6.2:  Let  us  again  consider  the  design  of  a  cantilever  beam.  Assume  the  beam  is 
subjected  to  a  point  load  w  at  the  free  end.  The  design  criteria  specifies  that  the  displacement  d  at 
the  free  end  does  not  exceed  the  a  specific  value  denoted  by  z.  From  mechanics  of  materials,  we 
know  the  performance  function  describing  the  end  displacement  is 


d  = 


wU 

3EI 


(6.77) 


where  w  is  the  load,  L  is  the  length  of  the  beam,  E  is  Young’s  modulus  of  the  beam  material,  and  /  is 
the  cross  sectional  moment  of  the  beam.  The  limit  state  that  separates  the  safe  region  from  the 
failure  region  is 


g{x)  =  z-d  =  2-  —  =  0  (6.78) 

^  ^  3EI 

Assume  the  deflection  criteria  z  is  deterministic  and  set  at  0.25  in.  Assume  the  load,  length,  modulus, 
and  moment  are  all  random  variables  with  the  following  distributions; 
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w~iV(/i^,cr^)  =  A^(5000,1000)  lbs. 
Z~iV(//i,<Tj  =  A^(l20,10)  in. 

E  ~  N{me,<^e)  =  A^(30X10^0.3X10‘‘)  psi 

/~A^{^,,o-/)  =  Af(l200,150)  in'* 

The  mean  deflection  according  to  (6.77)  is 

j  wL\  (5000)(120)^ 

3EI  3(30X10*)(1200) 


(6.79) 


Following  the  steps  outlined  above: 

STEP  1 :  Assume  the  MPP  is  at  the  mean  value. 


STEP  2:  Derive  each  of  the  derivatives 


X  =  (5000,1 20,30X1 0®  ,1 200) 


dx] 


Considering  the  random  variable  w ,  we  can  use  (6.51 )  to  write  (6.78)  in  terms  of  w'  as 

and 


3EI 


3EI 


Considering  the  random  variable  Z ,  we  can  use  (6.51)  to  write  (6.78)  in  terms  of  L'  as 


g  =  z-- 


and 


3EI 


^  -w3a^{aEL' 

dL'  3EI  El 


Considering  the  random  variable  E ,  we  can  use  (6.51)  to  write  (6.78)  in  terms  of  E'  as 

wL^ 

3{<7eE'  +  MeV 

and 

-wZ^(-lcr^)  _(JewL^ 

<3E'  3[aEE'  +  fiEfl  3Z'^/ 

Considering  the  random  variable  I ,  we  can  use  (6.51)  to  write  (6.78)  in  terms  of  /'  as 

_  wl} 

^~"~3E{ct,I’^Hi) 

and 
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_  -M'i^{-lcr,)  a/wL^ 

3E{(T,r+M,f  ~  3^:/' 


STEP  3:  Evaluate  each  of  the  derivatives  at  the  assumed  coordinates  of  the  MPP. 


dw' 

dL' 

dE' 


-(1000)(120y 

3(30X10*)(1200) 
-(5000)(10)(120)^ 
(30X10*)(1200) 
(0.3X1 0^)(5000)(1 20)^ 
3(30X10*)^(1200) 


=  -0.016 


-0.02 


=  0.0008 


^  (150)(5000)(120)^ 

dl'  3(30X1 0'^Xl  200)^ 

STEP  4:  Evaluate  each  of  the  direction  cosines  a*  according  to  (6.75). 


-0.016 


= 


-0.016 


V0.016^  +0.02^  +  0.0008^  +0.01^  0-0275 


=  -0.581668 


.  -0.02 

a,  - - 

^  0.0275 

.  0.0008 

.  0.01 

a,  = - 

‘  0.0275 


=  -0.727085 


=  0.029083 
=  0.363543 


STEP  5:  Write  each  of  the  coordinates  of  the  MPP  x*  in  terms  of  p  according  to  (6.76). 
w*  =  5000  +  (0.582)(1 000)yS  =  5000  +  582/7 
Z,*  =  120  +  (0.727)(10)/?=  120  +  7.27/7 
E'  =  30X10®  -  (0.029 l)(0.3X10®)y5  =  30X10®  -  873Oy0 
/*=  1200  -  (0364)(1 50)y9  =  1200  -  54.4y5 


STEP  6:  Substitute  the  x]  ’s  found  above  into  the  limit  state  equation  g(jt*)  =  0  such  that  the  limit 
state  equation  is  in  terms  of  p . 


Z{x\  025  (5000  +  582;g)(120  +  7.27/?)^ 

’  ■  3(30X10® -873O>9)(12OO-54.4y0) 


(6.80) 


STEP  7:  Solve  for  p  in  the  above  equation. 


Equation  (6.80)  has  multiple  solutions  for  /7 .  However,  we  are  usually  only  interested  in  values 

between,  say,  5  and  -5.  A  p  value  between  5  and  -5  corresponds  to  a  probability  of  failure  of 
between  0.0000003  and  0.9999997,  which  usually  covers  our  range  of  interest.  We  can  further 
'  reduce  our  range  of  interest  for  the  beam  example  by  noticing  that  the  mean  value  of  deflection  is 
0.08  inches  and  the  failure  criteria  is  0.25  inches.  More  often  than  not,  the  beam  deflection  will  be 
less  than  0.25  inches.  We  can  safely  say  that  the  probability  of  failure  is  less  than  50%.  Therefore, 
we  are  only  interested  in  positive  values  of  p .  Thus,  we  can  restrict  our  search  to  p  values  between 
0  and  5. 
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We  can  find  the  solution  to  (6.80)  by  trial  and  error  or  by  plotting  g{x)  vs.  p  as  shown  in  Figure 
26.  In  either  case,  we  find  that  p  =  3.66 . 


Figure  26:  g(x*)  as  a  function  of  p. 


STEP  8:  Find  the  updated  assumed  coordinates  of  the  MPP  x]  by  reevaluating  (6.76)  with  the  cr,  ’s 
and  p  found  above. 


w*  =  5000  +  (0.582)(  1 000)(3.66)  =  7130 

L'  =  120  +  (.727)(10)(3.66)  =  146.6 

£*  =30X10^  -(0.0291X0.3X10‘^)(3.66)  =  29.97X10* 

/‘  =  1200-  (0.05 1 7)(1 50)(3.66)  =  1 000 

STEP  9:  Repeat  steps  3  through  8  until  convergence  is  obtained  on  p . 

A  Summary  of  the  calculations  needed  to  converge  on  p  are  shown  in 

Table  5.  The  value  for  d  shown  in  the  first  column  of 

Table  5  are  the  most  probable  failure  point  found  by  evaluating  (6.77)  using  the  assumed  MPP 
values. 

After  iteration  #  2,  the  coordinates  of  the  MPP  x-  are  updated  by  reevaluating  (6.76)  with  the  a*  ’s 
and  p  found  in  iteration  #  2.  Each  of  the  derivatives  is  evaluated  at  the  MPP  and  each  of  the 
direction  cosines  a*  's  is  reevaluated  according  to  (6.75). 

Figure  27  compares  the  FORM  results  with  Monte  Carlo  simulation.  The  beam  tip  displacements  are 
plotted  as  a  function  of  the  probability  of  occurrence  with  standard  normal  units  p .  Even  though  the 
response  function  is  nonlinear,  it  can  be  seen  that  the  FORM  results  compare  very  well  to  the  Monte 
Carlo  Simulation,  which  can  be  considered  the  exact  solution. 
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Table  5:  Summary  of  iteration  results. 


Initial  Values 

Assumed  MPP 

Derivative  at  MPP 

Otl 

P 

w  =  5000 

L=  120 

E  =  30E6 

I  =  1200 
d  =  0.08 

w  = -0.016 

L  =  -0.02 

E  =  0.0008 

1  =  0.01 

w  =  -0.582 

L  =  -0.727 

E  =  0.291 

1  =  0.364 

3.66 

Iteration  #1 

MPP 

Derivative  at  MPP 

Otj 

P 

w  =  7130 

L=  147 

E  =  29.97E6 

1  =  1000 
d  =  0.2497 

w  = -0.0350 

L  = -0.0511 

E  =  0.00250 

1  =  0.0375 

w  =  -0.484 

L  =  -0.705 

E  =  0.0345 

1  =  0.517 

3.595 

Iteration  #2 

MPP 

Derivative  at  MPP 

otl 

P 

w  =  6738 
L=145 

E  =  29.96E6 

1  =  921 
d  =  0.2499 

w  =  -0.0371 

L  =  -0.0516 

E  =  0.00250 

1  =  0.0407 

w  =  -0.491 

L  =  -0.683 

E  =  0.0331 

1  =  0.539 

3.595 

FINAL 

MPP 

Derivative  at  MPP 

ai 

w  =  6766 

L=  144 

E  =  29.96E6 

1  =  909 
d  =  0.25 

w  =  -0.0370 

L  = -0.0519 

E  =  0.00250 

1  =  0.0413 

w  =  -0.487 

L  =  -0.683 

E  =  0.0330 

1  =  0.543 
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Tip  Displacement  (in.) 

Figure  27:  Comparison  of  FORM  results  with  Monte  Carlo  simulation. 


6.2.3  Mean  value  methods 

Mean  value  methods  are  particularly  suitable  when  a  closed-form  expression  of  the  limit  state  is  not 
directly  available,  as  is  the  case  when  finite  element  analysis  is  used  to  model  the  structural  response.  The 
reliability  computation  algorithm  used  in  the  rotor  hub  demonstration  problem  constructs  a  linear 
approximation  to  the  limit  state  in  the  equivalent  standard  normal  space.  Mean  value  methods  works 
differently,  as  follows. 

1 .  A  first-order  Taylor  series  approximation  is  constructed  at  the  mean  values  of  the  random  variables — in 
the  original  space — based  on  perturbation-based  sensitivity  analysis  using  the  finite  element  analysis  or 
other  implicit  response  models.  The  two-parameter  scheme  is  used  to  transform  this  closed-form 
approximation  and  the  random  variables  to  the  equivalent  standard  normal  space.  The  Rackwitz-Fiessler 
[21]  iteration  formula  is  used  to  estimate  the  MPP,  and  a  first-order  estimate  of  the  failure  probability  is 
obtained  similar  to  FORM.  This  is  referred  to  as  the  mean  value  first-order  (MVFO)  estimate. 

2.  The  MVFO  estimate  is  refined  using  advanced  mean  value  (AMV)  analysis.  This  is  simply  a 
deterministic  analysis  of  the  system  at  the  MPP,  to  re-evaluate  the  g-function.  If  the  MVFO  estimate 
were  accurate,  the  g-ftinction  would  be  exactly  zero.  If  the  MVFO  estimate  were  not  accurate,  then  a 
different  value  would  be  obtained  for  the  g-flinction.  The  MPP  and  the  probability  estimate  identified  by 
the  MVFO  are  now  assumed  to  correspond  to  this  new  value  of  the  g-flinction. 

3.  The  above  two  steps  are  repeated  for  different  values  of  the  g-function.  This  results  in  the  construction  of 
the  CDF  of  the  g-fimction. 

Mean  value  methods  are  a  practical  alternative  to  FORM,  if  the  deterministic  analysis  of  the  system  is 
expensive.  By  simply  combining  the  information  from  the  MVFO  step  and  one  additional  deterministic 
analysis,  one  is  able  to  obtain  a  substantially  improved  estimate  of  the  failure  probability. 

There  are  optimization  algorithms,  such  as  BEGS,  which  combine  the  information  in  the  iteration  steps  to 
perform  an  accurate  search.  Such  algorithms  are  especially  useful  when  the  Rackwitz-Fiessler  algorithm  fails 
to  converge.  However,  the  programming  effort  and  the  memory  storage  requirements  in  these  methods,  plus 
the  marginal  improvement  in  accuracy  over  the  Rackwitz-Fiessler  method  have  hindered  their  widespread 
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application  in  structural  reliability  studies.  The  AMV,  with  its  minimal  one  step  combination,  is  a  practical 
approach  in  this  direction. 

The  AMV  estimate  can  be  further  improved  with  more  iteration.  In  the  AMV  +  method,  sensitivity 
analysis  is  performed  after  the  AMV  step,  and  steps  1  and  2  of  the  above  procedure  are  repeated.  In  a  similar 
manner,  further  iterations  can  be  done  to  produce  AMV  +  +,  AMV  +  +  +,  etc.,  estimates.  In  each  case,  note 
that  only  the  first  iteration  (MVFO)  and  the  latest  deterministic  analysis  are  combined.  It  is  important  to  note 
here  that  the  sensitivity  information  in  the  AMV  is  computed  only  in  the  first  iteration,  i.e.,  at  the  MVFO 
step,  which  is  at  the  mean  values  of  the  random  variables.  If  the  sensitivity  values  are  needed  at  the  MPP, 
then  full  sensitivity  analysis  needs  to  be  performed. 

The  AMV  method  approximates  the  g  function  at  the  very  beginning  (at  the  mean  values  of  the  random 
variables).  That  is,  the  first-order  derivatives  are  calculated  at  the  beginning  using  finite  difference  or 
perturbation  analysis  to  have  a  first-order  Taylor  series  approximation  as  shown  below 

+  (6.81) 

where  x  represents  the  expansion  point.  Notice  that  the  linear  approximation  of  Eq.  (6.8 1 )  does  not  involve 
mixed  terms,  and  the  effects  of  the  random  variables  on  the  g-fiinction  are  considered  separately.  This 
introduces  additional  approximations. 


Numerical  Example  6.3:  Assume  we  have  a  cantilever  bean  of  length  L,  loaded  at  the  free  end  with  a 
force  w.  We  wish  to  find  the  vertical  displacement  at  the  free  end  of  the  beam.  We  know  that  the  load 
w,  stiffness  E,  and  the  moment  of  inertia  /,  are  all  normally  distributed  random  variables,  with  mean 
and  standard  deviation  as  follows: 

(5000, 1000)  lbs 

W  Z,«Af(120,10)  in 


Figure  28;  Cantilever  beam  loaded 
at  the  free  end. 

Let  us  assume  ,  for  this  example,  we  could  not  determine  the  displacement  in  closed  form  and  had  to 
use  the  simple  finite  element  model  in  Figure  28  to  determine  the  end  displacement.  The  details  of 
each  step  involved  in  performing  a  Mean  Value  First  Order  (MVFO)  analysis  for  the  above  example 
are  enumerated  below. 


£:«iV(30xl0%0.3xl0®)psi 
/«iV(1200,150)  in' 

We  know  from  mechanics  that  the  vertical  displacement  d  of  the 
free  end  of  a  beam  is 


cl  = 


wU 


(6.82) 


MVFO  Analysis 

Step  1 .  The  first  step  in  a  MVFO  analysis  is  performing  a  perturbation  analysis  about  the  mean 
values  of  the  input  random  variables  (w,  L,  E,  I)  using  the  finite  element  model,  to  determine  the 
vertical  tip  displacement  d,  of  the  cantilever  beam.  The  perturbation  analysis  involved  changing  each 
variable  by  one-tenth  of  its  standard  deviation  about  its  mean  value  and  rerunning  the  finite  element 
model  to  create  the  array  of  five  finite  element  runs  as  shown  below.  The  results  of  the  perturbation 
analysis  are  shown  in  Table  6. 
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Table  6:  Perturbation  analysis 


w 

L 

E 

/ 

d(FEM) 

5000 

120 

3.00E-r07 

1200 

8.00E-02 

5100 

120 

3.00E-t-07 

1200 

8.16E-02 

5000 

121 

3.00E-r07 

1200 

8.20E-02 

5000 

120 

3.00E■^07 

1200 

7.99E-02 

5000 

120 

3.00E-r07 

1215 

7.90E-02 

StepiZ  Based  on  the  perturbation  analysis  of  Step  1 ,  a  first  order  Taylor  Series  expansion  is 
performed  about  the  mean  values  of  the  random  variables  using  equation  (6.81 ),  expressing  the 
vertical  tip  displacement  d,  as  a  function  of  the  input  variables.  Equation  (6.83)  represents  this  linear 
functional  relationship.  Comparing  equations  (6.82)  and  (6.83),  it  shouid  be  noted  that  equation 
(6.83)  is  only  an  approximate  relationship  based  on  finite  element  analysis  at  the  mean  value  of  the 
random  variables. 


=  -0-083 1  + 1.6  *  1 0  *  w-f  2.1  ♦  1 0  ’  I  -  2.66  *  1  O’’  £  -  6.58  *  1  O'*  /  (6.83) 

Ste&a  Having  defined  equation  (6.83),  the  next  step  is  to  define  the  limit  state  equation  for  first  order 
reliability  analysis.  The  limit  state  equation  is  defined  as; 

S  =  (6.84) 

where  is  the  "level  of  displacement"  for  which  the  reliability  analysis  is  performed.  In  other  words, 

the  reliability  analysis  will  find  the  probability  that  the  beam  displacement  is  greater  than  d^ . 

Based  on  a  range  of  values  of  do,  the  probability  that  g  <  0  was  determined.  Table  7  lists  the  range  of 
values  of  do  considered  along  with  the  corresponding  reliability  indices  (P)  and  the  most  probable 
values  (MPP)  of  the  input  random  variables  (w,  L,  E,  I)  for  which  the  relationship  g  <  0  was  satisfied. 
Figure  29  plots  do  vs.  p  for  the  MVFO  analysis  which  represents  the  cumulative  distribution  function 
(CDF)  of  the  vertical  tip  displacement  dof  the  beam.  Notice  that  the  MVFO  curve  is  linear.  This  is 
because  the  random  variables  are  normal  and  a  linear  approximation  to  the  beam  displacement 
response  (equation  (6.83))  was  used  in  the  analysis.  However,  we  have  no  reason  to  believe  that  the 
true  displacement  response  of  the  beam  is  linear.  If  we  wish  to  develop  a  better  approximation  to  the 
true  CDF,  we  perform  an  advanced  mean  value  analysis. 

Table  7:  MVFO  results 


rfo 

P 

w* 

L* 

E* 

I* 

dpEM 

0.00 

3.20E■^00 

3.19E-f03 

9.62E-^01 

3.00E-^07 

1.37E-^03 

0.023 

0.02 

2.49E+00 

3.59E-r03 

I.OIE-^02 

3.00E-^07 

1.33E-^03 

0.031 

0.04 

1.78E■^00 

4.00E■^03 

1.07E-^02 

3.00E-^07 

1.29E-^03 

0.042 

0.06 

1.07E-r00 

4.39E-I-03 

1.12E-^02 

3.00E-^07 

1.26E-^03 

0.055 

0.08 

3.60E-01 

4.80E-I-03 

1.17E-^02 

3.00E-^07 

1.22E-^03 

0.070 

0.10 

-3.50E-01 

5.20E-r03 

1.23E-^02 

3.00E-^07 

1.18E-^03 

0.091 

0.12 

-1.06E■^00 

5.60E+03 

1.28E-^02 

3.00E-^07 

1.14E-^03 

0.114 

0.14 

-1.77E■^00 

6.00E■^03 

1.33E-^02 

3.00E-^07 

I.IIE-^03 

0.142 

0.16 

-2.48E-r00 

6.41E-t-03 

1.38E-^02 

3.00E-^07 

1.07E-^03 

0.175 

0.18 

-3.19E-I-00 

6.81  £•►03 

1.44E-^02 

3.00E-^07 

1."03E-^03 

0.219 

0.20 

-3.90E■^00 

7.21  £•►03 

1.49E-^02 

3.00E-^07 

9.95E-^02 

0.266 
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Figure  29:  MVFO  and  AMV  analysis  results. 


AMV  Analysis 


AMV  is  a  continuation  of  the  MVFO  analysis  in  which  we  use  the  MPP  from  the  MVFO  analysis  and 
the  finite  element  model  to  further  refine  the  reliability  analysis. 


Step  4.  The  vertical  tip  displacement  based  on  the  performing  a  finite  element  analysis  using  the 
MPP  values  of  the  random  variables  (columns  3  through  6  of  Table  7)  is  evaluated.  These  values  are 
listed  in  the  last  column  of  Table  7.  Notice  that  the  value  of  the  tip  displacement  from  the  finite 
element  analysis  is  not  the  same  as  r/,, .  The  value  of  is  equivalent  to  the  value  obtained  by 
using  the  approximate  equation  (6.83)  at  the  MPP.  The  displacement  value  at  each  p  Is  shifted  from 
d^  to  dppf^  as  show  by  the  AMV  curve  in  Figure  29.  The  AMV  curve  is  the  CDF  that  represents  a 
refinement  to  the  MVFO  CDF  using  just  one  additional  finite  element  run  at  each  value  of  p. 


A  Monte  Carlo  simulation  was  performed  using  the  closed-form  response  (equation  (6.82))  for 
comparison  purposes.  Figure  30  compares  the  AMV  results  with  Monte  Carlo  simulation.  The 
comparison  shows  excellent  agreement  between  the  two  methods.  However,  we  must  remember 
that  the  purpose  behind  performing  mean  value  methods  is  because  we  do  not  have  closed-form 
solutions  available  and  have  to  use  computational  models  to  determine  the  response  of  interest. 
Therefore  an  AMV+  analysis  can  be  performed  to  verify  the  validity  of  the  AMV  analysis. 


Figure  30:  Comparison  of  AMV  results  to  Monte  Carlo. 
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AMV4-  Analysis 


AMV+  is  a  continuation  of  the  AMV  analysis  in  which  we  use  the  MPP  from  the  MVFO/AMV  analysis 
and  the  finite  element  model  to  further  refine  the  reliability  analysis. 

Step  5.  If  we  are  interested  in  further  refining  the  AMV  estimate  at  the  p  value  of -3.19  (refer  Table 
7),  for  example,  we  perform  another  perturbation  analysis  by  repeating  Steps  1  and  2.  The  only 
difference  this  time  is  the  perturbation  analysis  is  now  performed  about  the  MPP  values  of  the  input 
random  variables  at  /9—3.19  as  opposed  to  the  mean  values  in  the  MVFO  analysis.  The  results  of 
the  perturbation  analysis  are  shown  in  Table  8. 

Table  8:  AMV+  perturbation  analysis 


w* 

L* 

E* 

/* 

d 

6.81  E+03 

1.44E+02 

3.00E+07 

2.19E-01 

6.91  E+03 

1.44E+02 

3.00E+07 

1.03E+03 

2.22E-01 

6.81  E+03 

1.45E+02 

3.00E+07 

1.03E+03 

2.23E-01 

6.81  E+03 

1.44E+02 

3.00E+07 

1.03E+03 

2.18E-01 

6.81  E+03 

1.44E+02 

3.00E+07 

1.05E+03 

2.16E-01 

Step  6.  Based  on  the  above  perturbation  analysis  a  first  order  Taylor  series  expansion  (equation 
(6.81))  is  performed  about  the  MPP  as  shown  in  equation  (6.85). 

=  -0.2266  +  3.2 1  *  1 0-^  w  +  4.588  *  1 0'^  Z  -  7.28  *  1  O'’  £  -  2.1  *  1 0'^  1  (6.85) 


Step  7.  The  limit  state  equation  is  defined  as  follows.  dAMv  is  chosen  as  0.21 9  because  at  the  AMV 
step  for  =-3.1 9  the  d  value  from  Finite  Element  Analysis  was  0.219  (refer  Table  7) 


S  —  ^  AMV 


approiiMii-, 


=  0.2\9-d. 


“PProx^n-^ 


(6.86) 


Step  8.  Performing  a  first  order  reliability  analysis  with  equation  (6.86)  as  the  limit  state  yields  the 
following  results. 


P 

w* 

L* 

E* 

/* 

0.219 

-3.41  E+00 

6.72E+03 

1.44E+02 

3.00E+07 

Figure  31  compares  the  AMV+  result  with  the  AMV  results.  It  can  bee  seen  that  the  AMV+  result  compares 
favorably  with  the  AMV  results  and  thus  it  can  be  assumed  that  the  AMV  and  AMV+  results  are  valid  in  the 
regime  of  p=  -3. 

6.2.4  Direct  FORM 

If  we  repeat  steps  1  through  8  the  analysis  would  be  known  as  AMV-H-.  If  we  again  repeat  the  steps  the 
analysis  is  known  as  AMV-H-+  and  so  on.  If  we  perform  the  analysis  in  which  the  number  of  repetitions  is 
based  on  the  convergence  of  p  and/or  the  MPP,  the  analysis  is  known  as  direct  FORM.  Therefore  MVFO 
and  AMV  methods  are  subsets  of  the  direct  FORM  analysis. 
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Figure  3 1 :  Comparison  of  AMV  results  and  AMV+  results. 

6.3  Design  of  Experiments  Input  Matrix 

Response  surface  methods  require  an  analytical  relationship  that  describes  the  component  behavior 
subject  to  a  specific  failure  mechanism.  In  cases  where  explicit  relationships  describing  system  degradation 
or  loading  exist,  they  can  be  used  by  probabilistic  methods  to  predict  reliability  and  life.  However,  if 
computational  methods,  such  as  finite  element  methods,  are  employed  to  estimate  component  loading  or 
behavior,  the  functional  relationship  between  loading  and  design  variables  must  be  approximated  using 
closed-form  equations. 

Response  surface  methods  calculate  implicit  functions  describing  component  behavior.  In  order  to 
estimate  the  response  surface,  several  steps  are  involved,  namely: 

1 .  Identifying  the  significant  random  variables. 

2.  Establishing  the  degree  of  model fidelity  required  or  desired.  The  engineer  must  determine  if  a  first- 
order  linear  model  description  of  the  phenomena  is  sufficient  or  if  a  higher  order  model  is  required. 

3.  Determining  the  experimental  design.  The  experimental  design  must  be  constructed  based  on  the  model 
that  must  be  fitted,  the  number  of  experimental  conditions  that  can  be  tested,  time  constraints  and  budget 
constraints.  In  addition,  the  number  of  levels  for  each  variable  must  be  determined  (i.e.:  how  many 
different  values  of  each  variable  will  be  tested). 

4.  Executing  the  experimental  design.  The  computational  software  should  be  run  at  the  individual  variable 
settings  (experimental  runs)  determined  by  the  experimental  design,  with  the  variable  settings  and 
resulting  output  recorded. 

5.  Fitting  the  response  surface  model.  Multiple  regression  techniques  are  used  to  estimate  the  parameters  of 
the  approximating  polynomial.  Typically,  first  or  second-order  models  of  the  response  surface  are 
sufficient. 

The  selection  of  the  model  order  is  fairly  straightforward.  Unless  evidence  or  knowledge  exists  to  the 
contrary,  a  second-order  design  is  desirable  since  first-order  effects  are  included,  as  well  as  measures  of 
curvature.  If  a  second-order  response  surface  does  not  sufficiently  fit  the  computational  model,  then  higher- 
order  designs  can  be  implemented. 

A  critical  consideration  is  the  number  of  computational  model  conditions  that  must  be  executed  in  order 
to  establish  a  valid  response  surface  model  of  the  phenomena  of  interest.  Minimizing  the  number  of 
experimental  conditions  that  must  be  modeled  can  drastically  reduce  the  computational  time  required.  While 
fractional  factorial  designs,  such  as  resolutions  III,  IV,  and  IV  designs,  might  seem  to  offer  the  most 
economical  means  of  determining  the  response  surface,  they  are  not  the  most  efficient  designs  available. 
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Figure  32:  CCD  treatments  for  2k  design  (left)  and  axial  runs  (right). 

In  selecting  the  experimental  design  for  fitting  a  response  surface  the  desirable  features  are: 

1 .  A  reasonable  distribution  of  data  points  throughout  the  design  space. 

2.  Enables  experiments  to  be  performed  in  blocks  to  reduce  error  and  noise. 

3.  Generates  estimates  of  error  for  the  fitted  model. 

4.  Economizes  the  number  or  runs. 

5.  Does  not  require  a  large  number  of  levels  of  the  variables. 

6.  Enables  determination  of  model  adequacy,  especially  estimation  of  lack  of  fit  of  the  model. 

7.  Allows  designs  of  higher  order  to  be  built  upon  sequentially. 

Although  fractional  factorial  design  reduce  the  number  of  combination  treatments  required,  they  require 
repetition  of  the  treatments  in  order  to  obtain  error  estimates  for  the  fitted  relationships.  In  contrast, 
experimental  design  tailored  to  fitting  response  surfaces  can  reduce  the  required  number  of  treatments  and  the 
number  of  required  repetitions.  If  the  factors  under  investigation  are  to  be  explored  at  two  levels  then  central 
composite  designs  (CCD)  have  several  advantages.  First,  CCD  enables  the  adoption  of  sequential 
experimentation.  In  this  case  a  2k  design  of  resolution  V  is  used  to  fit  a  first-order  model,  which  exhibits  a 
statistically  significant  lack  of  fit.  The  experiment  then  progresses  with  additional  treatments  to  include  2k 
axial  runs.  A  graphical  depiction  of  the  treatments  run  in  the  CCD  case  is  shown  for  2  variables  in  Figure  32. 
In  this  case  the  initial  runs  are  conducted  at  factor  level  settings  of  ( 1 , 1 ),  ( 1 ,- 1 ),  (- 1 , 1 ),  and  (-1,-1)  for  the  two 
factors.  If  following  the  initial  runs  the  fitted  first-order  model  lacks  fit,  additional  runs  are  conducted  at  5 
additional  treatments,  specifically,  (0,0),  (0,V2),  (0,-V2),  (V2,0),  and  (-V2,0)  [see  Figure  32].  The  axial 
treatments  [(0,V2),  (0,-V2),  (V2,0),  (-^2,0)]  can  involve  a  single  replication,  but  numerous  trials  are  typically 
conducted  at  the  center  point  (0,0).  The  addition  of  center  point  runs  allows  for  the  estimation  of  curvature 
effects,  while  axial  points  enable  estimation  of  quadratic  terms.  Thus  the  resolution  V  fractional  design 
allows  estimation  of  the  linear  and  2-factor  interaction  terms,  and  is  variance  optimal  for  these  terms.  The 
axial  points  enable  estimation  of  quadratics  terms,  since  without  them  only  pure  second  order  terms  are 
estimated.  Center  run  points  provide  estimates  of  error  and  improve  quadratic  term  estimation. 

Where  models  are  to  be  fit  to  factors  that  have  three  levels,  the  CCD  approach  is  not  applicable.  In  cases 
of  factors  with  three  levels  are  investigated,  Box-Behnken  designs  (BBD)  [24]  are  efficient  methods  to  fit 
response  surfaces.  Combining  2k  factorial  incomplete  block  designs  forms  the  BBD.  Examining  the  Box- 
Behnken  design  for  three  variables  in  Figure  33,  it  can  be  seen  that  it  is  a  spherical  design,  with  no  points  at 
the  upper  or  lower  variable  limits  (i.e.:  no  points  at  the  vertices).  This  can  be  advantageous  if  cost  or  physical 
constraints  make  conducting  tests  at  extreme  variable  levels  infeasible. 

There  are  other  advantages  to  using  design  specifically  tailored  to  generating  response  surface  in  addition 
to  that  of  sequential  experimentation.  A  significant  consideration  in  design  selection  is  that  of  rotatability. 

Design  rotatability  exists  when  N  Var  y{x)  /  is  constant  for  all  locations  equidistant  from  the  center  of 
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Figure  33:  BBD  design  for  three  variables, 
designs  for  establishing  the  treatment  array  for  the  experiments. 


the  design.  The  impact  of  rotatability  is  that  for  any  two 
points,  x,  and  Xj  that  have  the  same  distance  from  the 
center,  their  associated  predicted  values,  >’(X| )  and 
y(x2 ) ,  have  the  same  variance.  Thus  if  the  design  is 
rotated  about  the  center  the  variance  of  the  predicted  y  is 
unchanged.  The  implication  of  rotatable  designs  is  that 
all  predicted  response  values  are  equally  uncertain  at 
any  given  radius  from  the  center.  If  non-rotatable 
designs  are  used,  the  resulting  response  surface  could 
have  a  bias  in  the  variance  depending  on  both  the 
distance  from  the  center  and  the  factor  levels.  The 
benefits  rotatability  are  significant,  favoring  the  use  of 
either  rotatable  central  composite  or  Box-Behnken 


7  Phase  I  Conclusion 

Phase  I  of  the  proposed  research  effort  has  successfully  demonstrated  that  the  probabilistic  analysis 
framework  can  be  integrated  with  current  state-of-the-art  composite  design  to  estimate  the  failure  probability 
due  to  various  critical  failure  modes.  In  addition,  the  same  framework  can  be  used  to  compute  the  sensitivity 
of  the  various  random  input  design  parameters  to  the  final  design  without  any  extra  computational  effort. 

In  particular,  a  composite  helicopter  rotor  hub  test  specimen,  subjected  to  cyclic  loading,  was  analyzed  to 
predict  delamination  onset  failure.  A  structural  model  of  the  rotor  hub  was  created  using  the  finite  element 
program  ANSYS.  The  output  from  the  structural  analysis  in  ANSYS  was  fed  into  the  Virtual  Crack  Closure 
Technique  (VCCT)  module  to  model  delamination  onset.  A  Reliability  Analysis  framework  including 
Response  Surface  Analysis  and  FORM  modules  was  developed  to  integrate  with  the  structural  analysis 
modules.  The  response  surface  analysis  was  used  to  develop  a  limit  state  equation  relating  the  primitive  input 
design  parameters  (random  variables  such  as  E,  P,  and  O)  to  strain  energy  release  rate  G„,«.  FORM  analysis 
estimated  the  probability  of  delamination  onset  and  furthermore,  determined  the  distribution  of  the  failure 
probability  over  a  range  of  cyclic  lives.  In  addition,  the  sensitivity  analysis  in  FORM  revealed  that  the 
external  loading  to  the  structure,  namely  Z’  and  <[)  are  the  most  critical  parameters  affecting  the  reliability  of 
the  rotor  hub.  Summarizing,  the  following  tasks  will  be  completed  at  the  end  of  Phase  I. 

■  Demonstrated  probabilistic  analysis  framework 

■  Demonstrated  composite  Finite  Element  modeling  (FEM)  technique 

■  Developed  framework  for  integrating  probabilistic  and  FEM  codes 

■  Applied  developed  method  to  a  practical  problem 

Figure  34  further  helps  in  illustrating  the  tasks  accomplished  in  the  Phase  I  effort  and  how  these  tasks  fit  into 
the  complete  Phase  II  effort.  The  solid  boxes  in  Figure  34  represent  the  tasks  that  were  successfully 
completed  at  the  end  of  Phase  I.  The  hatched  boxes  in  the  figure  represent  the  tasks  that  are  proposed  in 
Phase  II  to  develop  a  powerful,  commercially  viable  composite  analysis  design  tool. 
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Figure  34:  Phase  I  and  Phase  II  tasks. 
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